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chaptir 1 


INTRODUCTIOIT: 

In any scientific st idy, proper description of the 
physical or chemical phenomenon for is an important part of 
the work. Hiese phenomena can gen Jr ally he expressed through 
mathematical models sho-^.'ing the relationship he .;ween the 
dependent and the independent variables. One expects that 
these models should he unique and which nplies that only one 
model should he sufficient to describe a oarticular phenomenon. 
Unfortunately the experimenter observes, specially in compli- 
cated cases, that a phenomenon can he described by more than 
one model. This is mainly hecc.use the mechaniems, on which 
the models ane based, are selcom perfectly kne^fn and consequ— 
ently the experimenter faces the problem of identifying the 
most adequate model, 

Por example in the case of non-newtemian fluid flow on0 
can obtain shear stress vs^ shear rate data which when fitted • 
to the different available fluid models e.g, Bingham plastic, 
Rheopectic, Thixotropic etc,, could give apparently meaning-! 
ful results in each case. Under this circumstances one has t^ 

properly classify the fluid, rather one has to identify the i 

[ 

most adequate model. ' 
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The numher of dependent variables in the fluid flow 
problem, explained above, is only one and is called a single 
response system. In chemical kinetics where more than one 
component is present during the reaction and observations are 
made on all the components the number of dependant variables 
is more thai one. In such a case the system is called a 
multiresponse one. 

While studying heterogeneous chema cal reaction which 
is a multiresponse situation, where most of the mechanisms 
have not yet been properly identified, one can propose more 
than one mechanistic model to describe the reacting system. Foi 
example while using Hou^gen and Watson's method of approach th( 
main steps involved are (i) Adsorption of the components on thf 
surface (ii) Reaction of the adsorbed components on the surface 
(iii) Desorption of the components from the surface. 

Considering any of the above step as the rate controllij 
step one can propose more than one mechanistic model for the 
same chemical reaction. These examples and many other which 
have not been explained necessitate* an ellaborate study of th 
methods for model discrimination. Also the use of multirespon 

I 

techniques for parameter estimation and model discrimination | 
not only increases the precision of the estimates of parameter 
but also reduces the volume of confidence ellipsoid. 
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The strategy of model discrimination ahd parajaeter 
estimation can "be described schematically as given below. 


I Computer) | Data 

__ 

An experimenter does his experiment in the laboratory 
which gives information in the form of data. These data aye 
used for payameter estimation and discrimination among rival 
models throxigh computer. The iterative process can be carried 
on until a satisfactory level is reached. 

LITERATURE SURVEY : 

Extensive literature is available on parameter estimatic 
and model discrimination for single response models but for 
multiresponse situations the technique available for both para-* 
meter estimation and model discrimination seem to be rather 
limited \jhich are discussed briefly in the following sections. 

Ray^^\ Singh and Rao^^^ have reviewed the techniques 

available for paiameter estimation in multiresponse situations, 
f 5) 

Box and Draper^ discussed the Bayesian estimation of para- 
meters. Draper and Hunter^ dealt with design of experimeh 
for parameter estimation in multiresponse cases. Beauchamp 
and Cornell^ mathematically converted multiresponse model 


I Labor atoryj^ 
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to single response one and gave a novel teclinique for simul- 
taneous non linear estimation of paraiaeters. Hunter^ 
reviewed different criteria for estimation of unknown constant 
from multire sponse data* Mezaki ap i Butt^®^ proved that 
Bayesian estimation of parameters is far better than genera- 
lized non linear least square technique. M.J. ^ ^ 

used the concept of non constant variance covariance matrix fo3 

( 12 ^ 

improved estimate of parameters. Box, et.al. ' discussed 
some problems associated with analysis of mtiltire sponse data* 

In model discrimination the same criterion used for 
discrimination in single response cases can be used for multi- 
response cases with proper modifications. These method 
comprise mainly of 

i) Likelihood methods 
ii) Bayesian methods 

Likelihood methods 

In these methods certain error structure is assumed a 
priori and based on this the likelihood of a particular model ; 
is Calculated. The likelihood ratios, defined as the ratios ' 

i 

of maxiinam likelihood of a given model to that of various 
other models are calculated and based on a likelihood ratio oi 
100 or more, discrimination is achieved. This method has bee 
amply exemplified by Rao, et.al.^^^^ 


in an earlier study. 
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Bayesian Methods 


In the case of Bayesian methods, certain prior proha- 
Mlitles axe assumed for the models snd in the light of the 
data the model probabilities are updated resulting in the post- 
erior probabilities. Invariably Bayes theorem is used for 
calculating the posterior probabilities for the different 
models as given by the relation 


pN+1 

r 



( 1 . 1 ) 


where is the posterior probability for the rth model after 

Jf 

(N+l ) runs and is the prior probability for the rth model 
before the (Bf+1 ) run. is the likelihood for the rth model an( 
V represents the total number of models considered. 


Also a knowledge of error structure is required for the 
calculation of likelihood and certain prior probabilities are 
to be assumed for each of the competing models. In the absence ; 
of any prior knowledge equal prior probabilities are assumed 
for all the competing models. Model discrimination is achieved ; 
when the posterior probability of any one of the models exceeds 
a preset value. 

Ray^^^ suggested a novel technique for model discriminati 

I 

by extending the parameter estimation criterion of Box and I 
Draper The method thcugh appealing failed to discriminate 


properly in some cases. 



6 


Sequential design for discrimination 

After conducting jj- experimental trials the posterior 
probahilit aes of the different competing models aie calculated 
as described in earlier section. If the discrimination is not 
achieved additional experimental ridis are necessary and the 
experiments are conducted at such conditions of the independent 
Variables which will give the best discrimination ajaongst the 
rival models* The criterion available for the sequential 
design of experiments are briefly discussed in the following 
sections : 

Roth’s^ criterion : 

For sequential design of experiments Roth suggested that 
the experiments be conducted at appropriate condition of 
independent variables which maximizes the objective function 
D given by 

“ = TT ( E IT ( V - P) 

i =1 r =1 1=1 

where i represents the responses, r represents the models and 
q represents the qth setting of experimental conditions. 
represents the prior probability of the rth model before the | 
(R+l ) experiment and ^^P^^esents the predicted value of the 

ith response for the 1th model and kth experimental setting. , 
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Box and Hill* criterion 

The most comprehensive technique that has so far been 
described is by Box and Hill, 

They suggested that the following objective function 
be maximized 



where all the terms and the underlying theory have been 
discussed in details in Chapter 2, 

( 1 7 ^ 

Hsiang and Reilly^ 'by suitable assumptions subverted 
Box and Hill’s assumptions of normality, known variance and 
linearity of models in the neighbourhood of parameter 
estimates* Prasad and Rao^^^^ for the first time used the 
expected likelihood instead of point likelihood in Box and 
Hill's criterion. They applied the modification to study 
catalyst fouling system and proved that the rate of discrimi- 
nation was much faster compared to the original Box and 
Hill's approach. 
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OBJECnVB OF THE PRESENT STUDY ; 

As seen from the a'bove literature survey the parameter 
estimation and model discrimination in mult ire sponse situations 
pose no speciai difficulties. However the effective use of 
sequeatial design of experiments has not been practically 
demonstrated barring for the only study by Prasad and 
In the present investigation it is proposed to see how 
effective is the 8eq[uential design of experiments for discri- 
mination, between rival models in multire^onse situations, 
proposed by Box and Hill^^^*^^^. Varioxis illustrative 
examples are considered to see how the posterior probabilities 
Vary during the course of the sequential design of 
experimentation. * 





chapter 2 

THEORY 
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2«1 Data G-eneratlon : 

In order to test tlie model discrimination criterion.^ ^ ^ 

for multiresponse case with unknown variance covariance 
matrix, with reference to Chemical Engineering problems, a 
chemical reaction was considered for the present study. 

A reacting system comprising of components A,B, and C 
undergoing a chemical reaction represents a general chemical 
reaction and some assumptions were imposed on the state of the 
reaction so that the study could be focussed on model discri- 
mination rather than on the detadls of the kinetic reactions. 

The assumptions are 

(1) The reaction is carried out in a batch reactor, 

(2) The reaction is isothermal* 

(5) There is no change in total moles, volume and density, 

(4) The reaction follows first order kinetics. 


The possible competitive models suggested are 


(i) 


(ii) 


(iii) 

(iv) 


A 


Q Q 

_ii^ B C 


6 


15 


A B C 

-! 2 ^ 0 




42 


^Q 


44 


( 2 . 1 ) 
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where 0 are the rate constajats or parameters and all the 
four assumptions made, hold a good for these models* 

The rate expressions for the models considered are 


Model 1 A 


B » C 


^11 =-dt^ = ®ii "i 


®11 ^ - «12 ^2 + ®13 ^3 


■^12 “ dt 111 12 2 

^ 13 “ " d ^ ®1 2 ^2 “ ® 1 5 ^3 


Model 2 


® 2 l ^23 

A r-^ B C 


21 “ dt 


- ©21 ^1 


’22 dt 


*21 “ *22 ■‘2 


*23 "^2 


d x, 

^23 ~ "“dt ®23 ^2 


®32 

Model 3 A ^ B -2f-> C 


d x^ 

^31 “TT 


'31 


^32 “ ~ S ^ ~ *31 " *32 ^2 

d x^ 

^33 “ “dTt ~ ^32 ^2 



Model 4 


®42 ®44 

Y4I = -dt = - ®41 ^42 ^2 

"^42 "dt^ ®41 ” ^42 ^2 ” ®43 ^2 ®44 ^3 

^43 " “ ®43 ^ ■ ®4I ^3 

where Y represents response, x represents concentration of 
components A,B,C and 0 represents rate constants or 
parameters. 

Due to the non availability of experimental data, at the 
time of this study, efforts were made to get simulated data 
conforming with the assumptions already made. The following 
represents the scheme adopted to help overcoming the above 
mentioned deficiency: 

(1 ) (a) Data were generated using model A-->B C 

(i) without correlation among responses, 

(ii) with correlation among responses. 

(B) Data were generated using model A B ^ C 

(i)without correlation among responses, 

(ii) with correlation among responses, 

(2) Data reported by (19) were considered. 

A comprehensive survey about generation of data has been 
( 1 ) 

reported by Ray but the data reported by Ray corresponds only 
to case 1B(i) and was also unrealistic in nature. In the 
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present study extensive data were generated corresponding 
to (1), compatible with the underlying theory given ly 
earlier workers^^^ representing real case study. 


Considering the general mvilti response model 

the data were generated after assigning some realistic values 


to 0. 


( 20 ) 


The method consists of calculating the x values from 


Thereafter f . was calculated from 

* w 


the analytical expression for x vs t for the particular model 
considered^^ ^ and at regular intervals of time *t‘ for the 
assumed set of parameter 
the rate expression and finally the normal error was added to 
give the simulated response voctor. To avoid the dependency 
in the expected values of responses'^ *' only (Iy2» ^r5^ 'wo.qre 
taken as the response vector as we have from material balahce 

+ B (1^2^ ® ^ ^ (2.1.1) 

Also €. , vector of errors for iiie qth experiment was 

w M. 

assumed to have normal distribution, i.e., 


E(e)=0 and £(€ , ^ ^ 

q — — q } 

(uxl ) fl xu) 


(uxu) 


( 2 . 1 . 2 ) 


Since the elements of € are jointl37 normal for simulation 

”^{ 1 ) 

Hie expression can be expressed' as a linear combination of 


independent random variables 
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' t, = . e (2,1.3) 

(uxl) (uxu) (uxl ) 

where is a non singular matrix 

such that = CT*^ C (2.1.4) 

and elements of e has a IT (0,1) distribution obtained from 
table of random normal deviates. While considering reaction 
kinetics study one inevitably faces the problem of choosing the 
response vector and the independent or controlled variable 
vector, j&n experimenter actually observes the concentrations 
of the components at regular intervals of time and it seems 
that it is more realistic to consider (x^, X 2 » ^ as the 

response vector and (x^ , x|, t)"^ as the controlled 

vector where the superscript o indicates initial state and 'f 
represents time. Under this circumstance it may be justified 
to assume that the error of observation will follow a multi- 
variate normal distribution. But it has been pointed out by 
( 1 

Ray^ '' that this consideration will not only make the system 
more complex but also may lead to noncovergence as far as 
Parameter estimation is concerned. 


( 


dt 


The other alternative is to consider the rate vector 
d Xp d x- 

» “d^ ' response vector and concentration 
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vector t ^ 2 * controlled variable vector* Bat 

since rates are generally not experimentally measurable 

quantities rather these are calculated from concentration vs* 

time data» the error distribution may not be normally 

distributed smA giving rise to non constant variance - 

( 9 ) 

covariance metrix discussed by Box^-^-', 

However, since simulated data for estimated value of 
concentration were available the second alternative was 
considered for data generation with the rate vector as response 
vector and also under thisc. circumstances it may be justifiably 
assumed that the error of observation will follow a multi- 
variate normal distribution. 

2*2 Estimation of parameters and variance - covariance matrix 
of observations 

The main problems excluding the inherent difficulties, 
associated with Box and Hill criterion'' ' for multiresponse 

situation are the estimation of parameters and prior knowledge 
of variance covariance matrix. 

( 22 ") 

Koch ' has discussed the estimation of variance in 

univariate case and also of covariance in multivariate case. 

( 23) 

E-ao' ' also discussed the estimation of variance covariance 

( 24 ) 

matrix from known data for responses. Hill and Hunter^ * 
extended Box and Hill criterion for single response systems. 
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analytically, for unknown variance case but no such method 

f •5') 

was available for mltire sponse case. Box and Draper^ 
suggested that covariance matrix may be replaced by cofactors 
of the elements in the dispersion matrix. But since the 
problem was toestimate parameters as well as covariaice matrix 
simultaneously, the method adopted by Beuchamp and Cornell^ 
for parameter estimation was used in the study. The inter- 
mediate covariance matrix generated while estimating parameters 
was used as the estimate for covariance matrix and the 
underlying theory is given as follows; 

Considering the multire sponse model 





If 

CA' V t> 

- = ^ ^ 1' 

(Uuxi ) 


/ 



fjH >’■ “'I 



( 2 . 2 . 1 ) 


( 2 . 2 . 2 ) 


are defined and it is assumed that for each j the n components 
of are generated from a distribution such that 

0 and (^^ I ^ ^ ^ (2.2.3) 

(Nxl) (Nxl) (Nx 1 )( 1 xU) (NxN) 

and also that the set of u random vectors follows a 

multivariate distribution with 

j - i ^ i where >0 and ij^j (2,2.4) 

(irx 1 )( 1 xiT) (ixN) 
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which implies that correlation between observations for 
different responses is being considered. 

There fere, 

E( ^ ) ~ ■-/I = 2- {'x' j (2*2, 5) 

(Fuxi ) (1 x3Ju) (NuxFu) (iToi) (#xlT) 

where X = -I; is a symmetric positive definite matrix and & 

r: 1 3 

denotes the Kronecker product, 

(Nxl ) 

are represented as shown and also it is assumed that ^ follows 
a multivariate normal distribution then equation (2»2,1) can 
be represented as 

Ij = (x,0)+ (2,2.6) 

Further, if T = (Y^ )‘^ (2.2.7) 

(full) ^ 2 ^ 

and F _ (F.(x,e)...P (x,©))"^ (2.2.8) 

(^uxl ) - ^ ^ 

then eqn. (2.2,1) can be more elegantly described as 

Y = F + 1 (2,2.9) 

(ijyXl) (NhXI) (Bill) 

which essentially resembles a single response model. 

With the above assumption the maximum likelihood estimate 


of parameters can be obtained by minimizing 
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-1 


(S.) = (I - (I - I) 


( 2 . 2 . 10 ) 


with respect to £, where eqn, (2.2,»10) is a function of £ 
only. The suggested estimate for 1. and hence is given as 


T (o) i i 

(12N) (NX1) 


(o) /N 


where ij(o) = (^^^(o) ••• 

(NX1 ) 

and ij|o) = , 4,(3)) 


( 2 . 2.11 ) 


where £ ^^^ is the least square estimate of £^^^ and is the 


vector of parameters for rth model and jth response and £ 

M 'i 

is a subset of £. The estimate for 0^'^'' is obtained ly 
non-linear regression of function by least square technique. 


(d) 


(e‘^b=(]Cj-F3{4j4))T(Yj-P.(5^.4) 


( 2 * 2.1 2 ) 


'3 ' '-3 

Modified Gauss-Newton algorithm was suggested for 
fitting of non linear regression and is ellaborated as follows 


P =1' + 8 n 

(hux ) (iruxl ) (Nuxm) 


(2.2.13) 


where jP' 

(Nuxm) 


(|^ • P^* ... 


(2.2.14) 


= {VV ' ^ 1 $ k $ m and f =ni!&i 2 ] 

Jx 


(2.2.15) 
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and 5 = £ - 0^5 (2.2.16) 

So that the least sq.uare estimate of $ is 

i = ((!*<; I) (2.2.17) 

( mxl ) (mxN'u) (MuxNu ) (Nuxm) (miNu) (UuxNu) (Nuxi ) 


Kie current estimate of vector © 

— o 

c © = © + c ^ -where o <r c < 1 

-0-0 - ^ -c 


(2.2.18) 


•when substituted in eq.n. (2,2.18) for obtaining the next current 
Value of O given by 


c ~1 


= © + c . S> 
-o mxn 


(2,2.19) 


M « 

and the iterative process is carried out until © is approximated 
with sufficient accuxacy. However, as other non linear 
regression methods^a feasible preliminary estimate of ©^ is 
necessary for convergence. 

In the present study all the models are linear with 

( 1 ) 

respect to parameters and method followed by Ray^ ' is 
adopted for estimation of parameters. 

After <T. . is obtained which gives the variance 
^ 3 

covariance matrix of observation between responses, -3 is 
obtained from eqn. (2.2.20) 

-5“ = 5 ® I (2.2.20) 

(nuxnu) (uxu) (nxn) 
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and estimated values of parameters given "by 

9= (X^ ^ 2)“^ )~^ X (2.2.21) 

This me'^hod-”. enaTales one to get simultaneous estimates of 

A , 

parameters 9 and i. which are essential in applying Box and 
Hill criterion for model discrimination. 


2 . ^quontial Design of experiments 


A general multiresponse situation represents the 
condition in which u multivariate observations or responses 
are being measured where the responses are functions of some 
independent variables, intrinsic parameters and, in addition, 
are also associated with the inherent error of observations. 


The above idea can be expressed mathematically in a 
more elegant way as follows 
For any model r, 




y^. = (2.3.1) 

where y*^ . is the jtb response for the q^th experimental condition, 

1!* J 

f^j(^j, 9y) is the expected value of y^^, 

is the vector of n independent variables for the 
qth experiment, 

0^ is the vector of the intrinsic parameter having 
m components, 

is the random experimental error of observation 
V St 

associated with the Jth response and the qth experiment. 
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The limits for the suffixes are 

(i) 1 ^r<v , where v number of models ore considered 

(ii) 1^3 ^u , where u n mher of responses are observed 

(iii) N ^ q. ^ H+H , where N set of observations haf4 already 

been made and design conditions for farther 
* 

" set are to be selected. 

N=0 implies that some prior information is 
already available. 

llso that 

^ ~ I where the parameters may or may 

not be same for all the models. 

(v) X = j 3^} ,l^.h < n 

(vi) The random error ^ . is assumed to follow multivariate 
normal distribution, and 
B(^2g^)=0 for all 

B(^iq_ for all i,3, 1=8 (2.5.2) 

= O for all i,3, q.5^s 

The last condition means that any two esperiments qL,s 
are considered to be independent of each other i.o., they are 
uncorrelated observations. 
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In vector notations, let 


Y — (y^ y^ ••• y^^) 


be the vector of observation, 

be the vector of the expected values 
of observation for the rth model, 


f 


T 




be the vector of the expected values 
of observation for all the models. 


It is assumed that the observations Y are normally 
distributed about their expected values f^ with a variance 
covariance matrix 

^ - -j ’}- f ^ Si - then the probability density for 


IB 


V JIxlL 



exp [- ^ f ) 


r' 

( 2 . 3 . 5 ) 


Now, if f^ is linear in the parameter space ©, or can 
be approximately expressed as a linear f -unction of 0 in the 
region of the parameter estimates 0, then one gets 


ra 






m 


z_ 

k= 


(©^ - 
'■^rk 


«rk)[ 




£2 




© 


rk 




«r='< 

-r -r 


( 2 . 3 . 4 ) 

»» 

■Where the ^ are estimated values of © based on the 
available If experiments (if N=0, then estimates of © must 
be known) • 



22 


m 


k=1 


where 






'rk 


9 =0 
*1— 


(2.3.5) 


( 2 . 3 . 6 ) 


and X 


^j.i 


j®. , 
r3,t 


y^ 

^r3,2 •■' ^rj,t!v. 


r3,2 


r3,tft- 


(2.3.7) 


Draper and Hanter^^’^^ have shown that the posterior 


density function for 0, after ]J statistically independent runs, 
is, 

P(© tly, 7) == L“ S)^M(0 - ^)J (2.3.8) 


where 


■u u ^ m 


(2.3.9) 


“ = X 71 - -4j Jrl 

j=1 1=1 ■ 

where is an element ofZ^ . 

This implies that (£. - 0) is normally distributed about 
0. With oovariance matrix . 

Por convenience, let 



- ^+1 
^r1,1 

#+1 
^r1,2 ■ 

yU+l -1 

r1,m 


- 


- ‘ f 


^+1 

^ru,1 

yN+1 



ru,2 

ru,m 

A 


( 2 . 5 . 10 ) 
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Then it has "been proved^ that (6 - ©) is 

normally distributed about 0 with a covariance matrix of 




(2.3.11) 


Also, because of the linearization of the models in the 
parameter space, f^ is normally distributed about the predicted 
response with a covariance matrix and the 

probability density function is given by 

- 1/2 


It- N ' r I t 1,^ i ir-i-KI,„ H-fK-1 

P(fr 1^“ 717^^“ ^ t ) (W^ ) 


(f -Y 
r r 


^ N+l 


)] 


( 2 . 3 . 12 ) 


which finally reduces to(l6) 

N+1 j -1 /2 


( 2 rr) 


(jS+1 _ 


where 


H+1 


Ir + 


H+1 


(2.3.13) 


(2.3.14) 


The following two relations with equation (2.3.13) gives 

= trace 


trace I^, 

is a (uxu) identity matrix 
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In Kullback 
hypothesis 


( 21 ) 


= trace T/ (i/ )-^ 

criterion for discrimination between two 


3+1 


(r:s)= (f;6) + ii-i) 


( 2 . 5 . 15 ) 




Vj7) 


where I^'^'^rrs) =j-|' p^) 


N+1 


-oo --oO 


= 11 ^ 


1+1 


N+1 


'I - trace I„4-trace 


u 


s 


.?/t1 )I( j.^N.1 )-1 , -^1^1 .‘^N+1 )J 


(2.3.16) 

A similar expression is obtained for (s;r) where the 

indexes of equation (2*5.16) are only interchanged. Maximizing 
y^'^\r:s) in eqn. (2,5«15) w.r.t, independent variables one 
gets the next suitable point for further experimentation. In 
this approach no weight wac given for the prior probabilities 
and Box and Hill removed the deficiency by maximizing a 
scalar discriminant function given as 



P^^jfKlM) 1(1:2) ... 1(1 :v)" 


V' 

1 1(2:1 ) 1(2:2) ... I(2:v) 

. 

: 


|j[(v:l) I(v:2) ... I(v:v) 
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where P denotes the prior probabilities for the models concerned 
and also that elements on the main diagonal of the I(r:s) 
matrix is zero# 

The expression for reduces to 

P^Pg J(1,2)+.P^P^ J(1,3)+ J{1,4)+P2P3 ^(2,5)+.,. 

+ P^,^ P^ J(v.1,v) (2.3.17) 

with eqn. (2.S.16) 

By use of equation (2.,5.17)/the expression for tahes 
the form 

? r r )■’ - )■’ 

r=1 s=r+1 

(t/+’ - (2.3.18) 

Bqn, (2.3.18) is used as a design criterion for choosing 
the independent variables for the (B+l )th run so as to effect 
discrimination between the rival models. 

The steps involved for discriminations among rival 
models using the discussed method of approach are as follows: 

(1 ) Experiments are carried out for N runs either based 
on judgement or by some arbitrary or saboptiimal way. 

(2) Based on these observed values of responses and the 
models under consideration, the parajaetcrs are estimated using 
linear or non linear regression techniques. 
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The coTaxiance matrix is generally unknown and so 
appropriate estimates are to he made (discussed in section 2.2), 

^■'5) Values of prior probabilities are to be calculated 
for (K+Oth run which are the posterior probabilities after 
the ITth run. To begin with the discrimination process, initial 
prior probabilities are assumed, based, on prior infomation and 
otherwise all are taken to be equal to 1/v. 

(4) is maximized and the corresponding set of 

independent variables is taken as the operating condition 

for the (N+l )th run, 

(5) Posterior probabilities after N4-1 run are calculated 
using suitable methods (discussed in section 2.5), 

(6) (i) If reaches some value such that it can be 

accepted by some criterion, experimentation is stopped, 

(ii) If the trend of change of is observed to 

follow some acceptable pattern, experimentation is stopped, 

(iii) If does not satisfy (i ) or (ii), experiment 

is conducted at ^+1 

and from step 2 onwards, all steps are 
are repeated until some model satisfies 6(i) or 6(ii), As the 
number of experiments goes on increasing one can delete models 
with very low posterior probabilities ahd also introduce other 
plausible models, depending on the judgement of the 
experimenter. 
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2.4 Maximization of the scalar disorlminanl; function : 

Optimization Technique 

In the present study the scalar discriminant function 
is to he maximized w.r.t. the vector of independent variables 
or controlled variables x. This essentially enables one to 
f in d the design condition x for the (1T+I)th run and ensixres 
that optimum discrimination takes place between the rival 
models . 

Efforts were made to standard optimization techniques 
but the following factors were considered before proper selection 
of the method was made. 

(A) Prior information of K . ; 

Since is a nultivariablo function, a knowledge of the 
nature of the response surface in the cu-ordinato space is 
essential because an initial point, x^, is necessary ■ 
for starting the search methods. Information of provides 
a feasible starting point and saves computer tine. No such 
prior information was available which could lead to an 
erroneous optimum point. 

(B) Derivatives of EL : 

The functional dependency of on the vector of 
independent variables x, was complicated in nature and no 
attempts were made to calculate the analytical derivatives of 
Ky. But attempts to calculate the numerical derivatives 
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using finite difference methods gave serious problems and the 
most important of them being the stability and selection of 
step size. Also, the advantage of the knowledge of numerically 
calculated derivatives are off-set by the time involved for 
computation* 

(C) Operability region : 

Though the maximization of a function gives the feeling 

of global maxima, in engineering problems, however, these are 

seldom achieved. This may be altributed to the fact that the 

Values of the independent variables predicted by optimization 

technique for a global maxima are not feasible in practical cases. 

For example it is not advisable to conduct a kinetic reaction at 
^and then to change the design condition to a temperature of -200°C 
a temperature of, say, 2000°CA,in the next run where both the 

design-conditions are predicted by theoretical considerations. 

In fact, while optimizing real problems, the co-ordinate 
space in which the function is to be maximized is bound by 
constraints, the values of which depend on the judgement of the 
experimenter. This space is called the operability region. 

Considering the factor (C) the following constraints 
were imposed on the independent variables 

(i) 0 :C x^ ^30 

(ii) 0 ^X 2 ^30 

(iii) 0 4 Xj ^ 30 


(2.4.1 ) 
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Where , X2 and x^ represent the concentration of AfB 
and C respectively in gm. moles/litre. These assumptions are 
justified in the case of kinetic reactions in liquid or gas 
phase, conducted in a laboratory scale batch reactor. 


Study of the function for different values of x were 
made only in the operability region, instead of the whole space. 
This study helped to choose a feasible starting point. However, 
factors (B) and (C) restricted the choice to get suitable 
derivative free multivariable constrained optimization technique. 

Complex algorithm by was used in the present 

study where the method maximizes a non linear multivariable 
function subject to non-linear inequality constraints. Mathe- 
matically, a function fCx^ , X2**- is maximised with respect 
to n independent variables x^ , X2»»».x^ subject to m constraints 
of the form 




(2.4.2) 


An initial point x® ic to be supplied and thereafter a 
complex of k^^n + 1 points is foimed inclusive of the initial 
point. The additional (k-1 ) points are generated using random 
numbers and the constraints for the independent variables are 
as 


^ij - % ^ij ^^i ■" Si) 


(2.4.3) 
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wiiBrs i. = 1# u 

3=1, .... k-1 , represents the vertex 

and r. - are random numbers distributed over the interval (0,1), 

1 0 

The points selected in this manner must satisfy the 
explicit constraints but may not satisfy the implicit 
constraints. If a point does not satisfy any of the implicit 
constraints, the point is moved half-way towards the centroid of 
the remaining points and this modification for the search 
point is continued until a satisfactory point is reached, 

x^^(new) = (old) + x^ ^JZ i=1,...l!r (2.4.4) 

where x^ "Zl i=1»...N (2,4.5) 

j=1 

The search technique is to compare the functional values 
over the constrained region or operability region. Initially 
the functional values are evaluated at each of the E vartices 
and the vertex which gives the lowest function value is 
discarded. The discarded point is replaced by another point 
times as far from the centroid of the remaining points 
as the reflection of the discarded point in the centroid. The 
new point thus obtained lines on the line joining the rejected 
point and the centroid, 

x^^(New) = - ^i j^°^^?)'^^ic' (2.4.6) 
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However, if the new point repeats itself in giving the 
least fimctional value then it is moved halfway towards the 
centroid to give a new trial point and is repeated until some 
constraint is violated. If che trial point violates any explicit 
constraint then the variable is reset to a value ^ = 0.001 inside 
the limit* Choice of and it depends on some factors which 
has "been extensively dealt with in the orginal work by 
and depends mainly on the experience and the nature of the 
function subjected to the appropriate constraints. The use of 
over reflection factor, o<>1, tends to cause an enlargement of 
the complex and thereby compensating for the shrinkage caused 
by the movement of the complex towards the centroid. Also if 
the initial point is far removed from the optimm it enhances 
the progress of the complex towards the optimum. Similarly 
k> n+1 helps the complex from collapsing and with simultaneous 
use of «^ > 1 helps the complex in maintaining its progress 
while it reaches some corner. 

The stopping criterion or the attainment of maxima was 
assumed when the function values of each vertex are within 
units for ^ consecutive operations. Also there was no 
systematic search for alternative maxima, the maximization 
was ensured by using more than one satisfactory initial 
point , 
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It may "be empliasized that complex method Toy Box 
provided a good search technique to get the constrained 
optimization, though not a perfect method. 

2 • 5 Posterior proba'bilities for the models 


The probabilities of all the competing models after the 
(H+l )th run was' updated using Bayes Theorem which is given below 





( 2 . 5.1 ) 


where is the posterior probability and is the prior 

probability of the rth model, being the likelihood of the rth 
model, V is the total number of competing models. 


If the values of the posterior probabilities of the models 
were below the accepted values the method of locating subsequent 
design condition was carried out via the steps for estimation of 
Parajneters and variance covariance matrix of observations 
until an acceptable value of the posterior probability for a 
Particular model was attained. 
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CHAPTER 3 

RESULTS MB DISCUSS IONS : 

5*1 Data Generation . 

Nine sets of data -were considered to test the Box and Hill *s 
criterion for model discrimination in miltiresponse situation. 

The scheme adopted for data generation is given helow; 


( 1 ) 


( 2 ) 


Model Used s C 

Data Set 1 


. -1 . -1 

(i) I^rameter values : ~ 

(ii) C matrix (eqn 


.2.1.3)s /. 005.0 

\.o .005 J 


Data Set 2 


(i) Ifeirameter values 

(ii ) C matrix 

sp 

Data Set 3 

(i) Parameter values 


e_. = 0.1 min”"' = 0.1 min 

51 52 

/ .003 .005 

( ,003 .003 


=0.2 min ^ ©,„ = 0.1 min""' 

3 ' 32 


(ii) C 

matrix 

• 

0 

0 

VM 

.0 




0 

0 

Data Set 4 



(i) Ihrameter values * ©^^ 

= 0.2 min 

-1 

(ii) C 

matrix : 

/ .005 

.003 



{ .005 

.003 


©. 


32 


0.1 min 


Data Set 5 

Data reported (19 ) using model A — >B 
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(3) mel Used : A®11— 


6 


15 


Data Set 6 s C matrix 

I^rameter values : 

(ii) C matrix 

Data Set 7 

(i) Parameter values t 

(ii) 2 matrix : 

Data Set 8 . 

(i) I^^rameter values t 

(ii) C matrix t 

Data Set 9 . 

(i) Rarameter values : 

(ii) g matrix : 

luring data generation, 
maintain, the jbysical behaviour as far as possible. !Ehat is corresponding to 
an increase in concentration of a particular component, the response udiich 
is the rate of change of concentration should be positive and vice versa. 
However in practice this consideration doesndt hold good specially for small 
changes in concentration and hence for smaller values for rates. Thus to 
simulate real situation, it was observed that a value of O.OO5 to O.OO5 for 
the elements of C matrix sufficed the requirement. The logic was that the 
error introduced with the actual rate was of the form C e_, yAiexe e_ is the 
vector of random normal deviates and the largest element in e was roundabout 


j .005 .0 ^ 

v° - 9 °v .1 _1 

©^^=0.08 min 6.^2=. 02 min ©.^^=.06 min 
.005 .0 
.0 .005 


©.^^= .04 min ^ ^ 

/.005 .0 
l.o .005 J 


“1 -1 -1 
.02 min min ©.j^=,C1 min 


.005 .005 \ 
.005 .005 I 


©^.j= .05 min ©.^2« .01 min~^ ©^^=.01 min~^ 

.005 .005 
.005 .005 


the main consideration was that the data should 
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+ 2.5> and thus the magnitude of largest error that can occur was nearly 
+ .008 to + .013, It should be mentioned that in the present study the response 
generally Taried from .8 to .001 and thereby the higher rates were less affected 
tdiereas the lower rates were affected increasingly, justifying the assumptions 
made. 

Initially I5 runs were generated using the scheme as already explained 

and are shown in Tables '^.1 through Table 5*9 iiihere represent the 

concentration of A,B and C respectively in gm moles/litre and Rate 2, Rate 5 
the responses of B and G respectively in gm moles/litre min. 
represent/ Bie models given be equation 2.1 were considered as the competing 

models and the parameters and variance covariance matrix of observation were 

next estimated. 

5 *2 Estimation of Ihrameters and Variance -Covariance matrix. 

The method used for the estimation purpose was developed by Beauchamp 

and Comell^^^. The computer program used for the above purpose was developed 

( 1 ) 

by Ray' and proper modifications were made in it to suit the present study. 

The estimates for 0 and Xr are given in Tables ^.10 through 

As seen from Table 3 corresponding to model 3j data set 1 and run 

15 » the estimated values of the parameters are 0.1 OC4 min and 
. 1 

.0998 min which are very close to the actual values of the parameters 

.- 1.-1 

.1 min and .1 min , used for data generation. But for the estimates of 
^ r, the ^^22 element i.e. the 4th element in the column for hr is seen 



Table 3.1 
DATA SET 1 
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0 * 7 # 

Model = .1 min 0,^ = .1 mm C = (.005 0 ^ 

“ O .003^ 


ROT 

Ro. 

_ f 1- - 

^2 

- 

RATE 2 

RATE 3 

1 


6.703 

2.681 

0.615 

0.405 

0.271 

2 


4.493 

3.595 

1.912 

0.081 

0.357 

5 


3.012 

3.614 

3.374 

-0.056 

0,366 

4 


2.019 

3.230 

4.751 

-0.118 

0.326 

5 


1.353 

2.707 

5.940 

-0.134 

0.272 

6 


0.907 

2.177 

6.916 

-0.128 

0.217 

7 


0.608 

1.703 

7.689 

-0.113 

0.166 

8 


0.408 

1.304 

8.743 

-0.089 

0.131 

9 


0.273 

0.984 

8.288 

-0.071 

0.098 

10 


0.183 

0.733 

8.743 

-0.063 

0.066 

11 


0.123 

0.540 

9.084 

-0.036 

0.059 

12 


0 . 082 

0.395 

9.337 

-0.031 

0.040 

13 


0.055 

0.287 

9.522 

-0.019 

0.033 

14 


0.037 

0.207 

9.756 

. -0.015 

0.023 

15 


0.025 

0.149 

9.826 

-0,016 

0.011 

16 


0.001 

29.999 

0.001 

-2.991 

2.995 

17 


0.000 

29.999 

3. 099 

-3.006 

3.011 

18 


0.001 

29.999 

18.579 

-2.998 

3.005 

19 


0.001 

29.999 

18.579 

-3.006 

3.012 

20 


9.222 

0.001 

29.999 

0.936 

0.004 

* 

X 

are 

in gm moles/1 ,. 




Rates 

axe in 

gm mole /I. 

..min. 
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lable 3.2 


^31 

Model I A — > 

®32 

' B ’“^C 

MIA 

SET 2 

rQin~^ ®52 “ 

.1 0 

o o 
o o 

• • 

oo 
o o 

• • 

7 

RIM Mo . 

^1 

^2 

^3 

Rate 2 

Rate 3 

1 

6.703 

2.681 

0.615 

0.406 

0.272 

2 

4.493 

3.359 

1.912 

0.090 

0.359 

3 

3.011 

3.614 

3.373 

- 0.055 

0. 366 

4 

2.019 

3.230 

4.750 

- 0.123 

0.320 

5 

1.353 

2.706 

5.939 

- 0.130 

0.275 

6 

0.907 

2.177 

6.915 

- 0.125 

0.219 

7 

0.608 

1.702 

7,689 

- 0.108 

0.171 

8 

0.407 

1.304 

8.288 

- 0.090 

0.129 

9 

0.273 

0.983 

8.743 

- 0.070 

0.098 

10 

0.183 

0.732 

9.084 

- 0.060 

0.067 

11 

0.122 

0.540 

9.337 

- 0.035 

0.060 

12 

0.082 

0.395 

9.522 

- 0.028 

0.041 

13 

0.055 

0.286 

9.658 

- 0.018 

0.033 

14 

0.037 

0.207 

9.755 

- 0.015 

0.022 

15 

0.024 

0.148 

9.826 

- 0.017 

0.010 

16 

0.001 

29.999 

3.799 

- 3.031 

3.035 

17 

0.001 

29.999 

17.382 

- 3.060 

3.065 

18 

0.001 

29.999 

29.657 

- 3.045 

3.051 

19 

0.001 

29.999 

0.001 

- 3.060 

3.066 

20 

25.713 

0.001 

29.999 

2.690 

0.004 
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Table 3.3 
DATA SET 3 


([odel : 

®31 S2 

A— ^ B — > 0 

031 - 

.2 inirL~^ ®32 “ 

.1 S = (‘“o^ 

iun Ho. 

^1 

^2 

5 

Bate 2 

Bate 3 

1 

4.493 

4.420 

1.087 

0.459 

0.444 

2 

2.019 

4.949 

3.032 

-0.093 

0.493 

3 

0.907 

4.209 

4.883 

-0.235 

0.425 

4 

0.408 

3.223 

6.370 

-0.238 

0.325 

5 

0.183 

2,340 

7.476 

-0.196 

0.235 

6 

0.082 

1.105 

8.208 

-0.150 

0.164 

7 

0.037 

1.142 

8.821 

-0.111 

0.110 

8 

0.017 

0.782 

9.201 

-0.075 

0.078 

9 

0.007 

0.531 

9.461 

-0.051 

0.053 

10 

0.003 

0,360 

9.657 

-0.043 

0.028 

11 

0.001 

0.242 

9.756 

-0.019 

0.030 

12 

0.001 

0.163 

9.836 

-0.016 

0.017 

13 

0.001 

0.110 

9.890 

-0.007 

0.015 

14 

0.001 

0.074 

9.926 

-0.005 

0.009 

15 

0.001 

0.049 

9.950 

-0.009 

0.001 

16 

0.001 

29.999 

0.001 

-3.021 

3.026 

17 

0.001 

29.999 

29.999 

-3.046 

3.051 

18 

0.001 

29.999 

3.470 

-3.033 

3.039 

19 

0.001 

29.999 

29.999 

-3.046 

3.053 

20 

29.999 

7.132 

29.999 

-5.592 

0.726 
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Table 3.4 


Si S2 

Model : 9 


31 


DATiA SET 4 


2 min~^ 0^2 = •! min' 


Run Eo. 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 
11 
12 

13 

14 

15 

16 

17 

18 


4.493 

2.019 

0.907 

0.407 

0.183 

0.082 

0.037 

0.016 

0.007 

0.003 

0.001 

0.000 

0.000 

0.000 

0.000 

29.999 

0.001 

0.001 


4.419 

4.948 

4.209 

3.222 

2.340 

1.649 

1.142 

0.782 

0.531 

0,359 

0.242 

0.163 

0.109 

, 0.073 
0.049 
0.001 
29.999 
29.999 


X5 


1.086 

3.032 

4.883 

6.369 

7.476 

8.267 

8.820 

9.201 

9.461 

9.637 

9.756 

9.836 

9.890 

9.926 

9.950 

29.999 

10.577 

13.073 


Rate i. 


0.460 

-0.091 

-0.234 

-0.243 

-0.192 

-0.146 

-0.105 

-0.075 

-0.051 

-0.041 

-0.017 

-0.013 

- 0.006 

-0.005 

-0.009 

6.006 

-3.004 

-3.004 


n _ /.003 .003^ 
~ " ^003 .003^ 


Rate . 


0.446 

0.495 

0.425 

0.319 

0.238 

0.166 

0.115 

0.077 

0.053 

0.030 

0.030 

0.018 

0.015 

0.009 

0.000 

3.049 

3.009 

3.010 
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Model ; 

(19) 

ABC 

Table 

DATA 

3.5 

SET 5- 



Run Ho, 

^..1 

X2 

^3 

Bate 2 

Rate 3 

1 

3.862 

0.137 

0.000 

1.372 

0.008 

2 

3.726 

0.271 

0.003 

1.337 

0.018 

3 

3,464 

0.525 

0.009 

1.276 

0.034 

4 

3.229 

0.744 

0.026 

1.093 

0.082 

5 

3.011 

0.941 

0.047 

0.983 

0.108 

6 

2.807 

1.118 

0.074 

0.884 

0.134 

7 

2.617 

1.277 

0.105 

0.793 

0.156 

8 

2.440 

1.419 

0.140 

0.710 

0.175 

9 

2.121 

1.557 

0.220 

0.596 

0.200 

10 

1.843 

1.845 

0.311 

0.466 

0.228 

11 

1.603 

1.985 

0.411 

0.358 

0.249 

12 

1.593 

2.089 

0.517 

0.259 

0.265 

13 

1.129 

2.810 

0.690 

0.152 

0.288 

14 

0.915 

2.208 

0.876 

0.046 

0.311 

15 

0.740 

2,211 

1.048 

0.005 

0.286 

16 

0.001 

29.999 

29.999 

-3.558 

3.563 

17 

0.009 

3.113 

29.308 

-3.849 

3.933 
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Table 5.6 


^11 %2 
Model : A B 

®13 


DATA SET 6 


11 

12 

i3 


.08 min 

r 0 • “1 

.02 mm 
.06 mil“^ 


n _ r.005 

= - Ko 


.0 

.005 


) 


Sun So, 

^1 • 

X2 

23 

Sate 2 

Bate 3 

1 

10.892 

5.295 

1.813 

0.875 

- 0.002 

2 

7. ,909 

8.174 

1.916 

0.579 

0.043 

3 

5.745 

10.092 

2,164 

0.385 

0.070 

4 

4.171 

11.359 

2.470 

0.256 

0.081 

5 

3.028 

12.188 

2.784 

0.161 

0.073 

6 

2.199 

12,723 

3.078 

0.106 

0.070 

7 

1.597 

13.064 

3.340 

0.064 

0.058 

8 

1.160 

13.276 

3.565 

0.045 

0.056 

9 

0.842 

13.405 

3.753 

0.027 

0.045 

10 

0.611 

13.480 

5.909 

0.006 

0,027 

11 

0.444 

13.521 

4.035 

0.013 

0.034 

12 

0.322 

13.541 

4.137 

0.004 

0.023 

13 

0.234 

13.548 

4.218 

0.004 

0.022 

14 

0.170 

13.549 

4.281 

0.002 

0.017 

15 

0.123 

13.545 

4.331 

0.004 

0,017 

16 

29.999 

0.000 

29.999 

4.315 

-1.913 

17 

0.001 

29.999 

29.999 

1 .210 

-1 .205 
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Sable 3*7 


DATA SET 7 


Model : 

®11 ®12 
A— ^ 

01 5 

Qll = 
®12 = 
013 = 

OJiX i 

—1 

,04 min 

T • ~1 

,1 min 
. 2 min’"^ 

p, _ /.005 

S: -k qt 

.JO X 

• *^ao5 ^ 

Eun Mo. 

^1 ' 

^2 


Rate 2 

Rate 3 

1 

12.782 

3.601 

1.617 

0.475 

0.037 

2 

10.892 

5.137 

1.971 

0.310 

0.114 

3 

9.282 

6.240 

2.478 

0.240 

0.125 

4 

7.909 

7.118 

2.972 

0.200 

0.119 

5 

6.740 

7.848 

3.412 

0.163 

0.098 

6 

5.743 

8.464 

3.792 

0,142 

0.088 

7 

4.894 

8.988 

4.118 

0.117 

0.072 

8 

4.171 

9.433 

4.396 

0.107 

0.068 

9 

3.554 

9.813 

4.633 

0.089 

0.057 

10 

3.028 

10.136 

4.835 

0.067 

0.038 

11 

2.581 

10.412 

5.007 

0.069 

0.045 

12 

2.199 

10.647 

5.154 

0.054 

0.034 

13 

1.874 

10.847 

5.279 

0.050 

0.032 

14 

1.597 

11.017 

5.386 

0.042 

0.027 

15 

1.361 

11.163 

5.477 

0.039 

0.026 

16 

29.979 

0.000 

30.000 

7 .311 

-6.110 

17 

0.001 

29.999 

0.001 

-3.011 

3.01 6 
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Table 3.8 
DATA SET 8 

®11 

Model : A-H. B*;:0 6 ^ = .02 0 = (;g05 . 005) 

3 9^2 ^ 2 mlrL 

8„ = .01 min-^ 


Run No, 

^1 

X2 

^3 

Rate 2 

Rate 3 

1 

7.686 

0.346 

1.967 

0.171 

-0.008 

2 

7.095 

0.962 

1.942 

0.145 

0.002 

3 

6.550 

1.487 

1.963 

0.121 

0.010 

4 

6.046 

1.933 

2.021 

0.105 

0.021 

5 

5.581 

2.310 

2.108 

0.084 

0.023 

6 

5.152 

2.628 

2.220 

0.072 

0.030 

7 

4.756 

2.894 

2.350 

0.062 

0.036 

8 

4.390 

3.116 

2.493 

0.052 

0.039 

9 

4.053 

3.300 

2.648 

0.042 

0.040 

10 

3.741 

3.450 

2.809 

0.034 

0.041 

11 

3.454 

3.572 

2.974 

0.029 

0.043 

12 

3.188 

3.670 

3.142 

0.023 

0.043 

13 

2.943 

3.748 

3,309 

0.018 

0.043 

14 

2.717 

3.807 

3.476 

0.012 

0.040 

15 

2.508 

3.852 

3.640 

0.005 

0.037 

16 

0.001 

29.999 

27.757 

-0.330 

0.335 

1 7 

29.999 

0.001 

29.999 

0.880 

-0.257 
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Table 3.9 


Model 


MTA SET 9 

®11 ®12 _ 
A — 0,-, = .03 min” 

n m • “1 

13 0^2 = .01 min 

= .01 min”^ 


n _ /•005 
- ■ ^.005 


.005x 

.005^ 


Run No, 

^1 

^2 

1 

9.418 

2.537 

2 

8.353 

3.485 

3 

7.408 

4.283 

4 

6.570 

4.952 

5 

5.827 

5.508 

6 

5.169 

5.968 

7 

4.584 

6.344 

8 

4.066 

6.650 

9 

3.606 

6.894 

10 

3.198 

7.086 

11 

2.837 

7.234 

12 

2.516 

7.345 

13 

2.231 

7.423 

14 

1.979 

7.476 

15 

1.755 

7.505 

16 

17.164 

11.897 

1 7 

0.001 

29.999 


^3 

Rate 2 

Rate 3 

0.045 

0.262 

0.029 

0.162 

0.220 

0.036 

0.308 

0.183 

0.040 

0.478 

0.155 

0.047 

0.665 

0.124 

0.046 

0.864 

0.103 

0.050 

1.072 

0# 086 

0.054 

1.285 

0.070 

0.056 

1.500 

0.055 

0.055 

1.716 

0.042 

0.053 

1.929 

0.034 

0.055 

2.139 

0.024 

0.053 

2.345 

0.017 

0.052 

2.549 

0.009 

0.048 

2.739 

0.001 

0.044 

23.085 

0.649 

-0.123 

29.999 

0.032 

-0.027 
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5 

to be nearly 10 times greater than the remaining elenents and also the 

, .Qlement 

estimated value for the (2,2}/wa8 not in good agreement with the value used 
for data generation. 

5.5 Sequential Design of Experiments. 

O- j. 

After obtaining the estimated values of © and Zr fhe scalar % 

( 25 ) 

discriminant function Kv was maximized using Box's algorithm. 

She elements of ( ir j formed corresponding to the highest value 

of Kv within the operability r^icaa are reported in. Tables 3-16 through 5*24. 

Also tihe values of the independent variables x^* x^ namely concentrations 

of component A,B and C are also presented in those tables. These values of 

"til 

the independent variables were used for ccnducting the (u+l ) run i.e. the 
next experiment and the values of the responses were - simlated using the same 
strategy as that of the original data generation. These concentration data 
with the respcsises for different sets are shown in Tables 5.I through 3.9* 

3.4 Posterior Erobability . 

After conducting the (N+l ) run the posterior probabilities of all 
the models were calculated using equation 2.5.I and equation 2.3.13* 

If the posterior probability was below .9 for the models the idiole 
sequence firom 5 *2 to 3*4 "was repeated and the probabilities are reported in 
Tables 5«16 through 3.24* 

An overall review of the Tables 3*16 through 3*24 indicates that 

(i) For data sets 1,2,3 the posterior probabilities for the correct 

models (i.e. Model 3) showed an increase from .25 to O .48 nearly after 
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0equential design. In the subsequent runs the probabilities increased 
gradually but wilh some oscillations to values round about 0.6. 

(ii) For data set 4» the posterior probability for the eojrrect model 
(model 5) increased from ^25 to .52 after the 16th run. However during 
the 19th run the probability reduced to 0.53(heai‘ly) end further 
calculation was stopped. Explanatiott is given in the next paragraph. 

(iii) For data sets 5>6,7»8»9 after the 16th run i.e. for the 1st sequential 
mm the posterior probabilities for the correct model showed a rapid 
increase from the initial model probability. That is from the initially 
assigned model probability of .25 it went as high as .8 to .9. After 
iTth run the probabilities of the correct model approached 1 and 
further calculation was not necessary. 

Bie expression that is used for finding the posterior probabilities 
of the models is given by equation 2.5.1 idiere the likelihood for the models 
after the run is calculated using equation 2,5.13. Ihe most dominant 


factor in equation 2.5*13 was the exponential factor and the exponent on this 
factor cmtains expression involving (ITr^"^^) and )• 5hrom 

the tables 5*16 through 3*24 it can be seen that the ( 2.^ j matrices are 


ill conditioned for most of the cases, as indicated by the small values of 


the determinant for . Under tiiis circumstance*, if the difference 

N+l H+l \ 

(Y ” Yr very small, the likelihood for the models approaches 

very small values. However ^in the present study the error added for the 
correct model was small in magnitude and didnot pose any serious problems for 


discriminatictti pm*pose. 
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Looking into the posterior prohahilitgr values for data set 1 , 2 and 3 » 
it oan he seen that model 1 and model 3 (model 5 Is the correct model ) are 
eoay?etiag Pereas the remaining model probabilities approaoha very small values. 
The main reason is that during parameter estimation the parameter for ■tJie 

•7 

1 at model approached zero and virtually there was no difference between 
model 1 and model 5 except for liie estimates of variance covariance mateix 
for obvitws reasons and consequently -fee difference between ) 

was much smaller for model 1 and 5 compared to model 2 and 4 • These factors 
aooount for the probability values as obtained which is expected. Therefore 
Box and Hill *8 criterion was adequate enough to discriminate between »ival 
models for data set 1,2 and 5* 

For data set 4 frm Table 3 *19 it is seen that the (1 ) element of 

( for the I9tii arun was much smaller for model 1 compared to model 5» 

Also the values for (Y - Y ) for both models are comparable with each 

T 

other as explained above, and hence model 1 (the incorrect model) prevailed 
over model 3 (the correct model). 

For data sets 6,7 »Q and 9 as already pointed out, the values of 
the model probabitiea ^owed a large change after the 1st sequential run. 

Bie reason bdbiind this is that the difference (i -Y^ ) was significant 

for all the models except model 1 (the correct model) and hence the likelihood 
values for model 1 most, have taken lai^ger magnitudes compared to other models 
accounting for the very high posterior probabilities as shown in the Tables. 

For data set 5, the correct model is 3 and the same analysis as 
given for data sets 6,7, 8,9 holds good in this case. 
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Here it should also be pointed cut that in the expression for 
equation 2.3*18 , ■feough the prior probabilities of the models are used 
as weights but as in ihe present case nhere the ( 27^. )*" has large 

elements, the weights attached can be offset by the remaining factors 
containing ( ) terms. This may lead to inefficient design and may 

cause fluctuations in the posterior probability values. 
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lllie probability density function of Y 


N+1 


j B+l 1 

1 V 1 l.^r 1 

-z 

r 

! ^ " {2nT 

- exp 


given -2.^ IS 


^T 


)-'> (y 1T+1_ 9^N+in...(4.l) 




Tidiere i r 


N+1 


r 


r + W 


N+1 


(4 *2) 


N+1 

If the marginal probability density function of Y^ is considered then 

00 

P (y/'^^I 2: P ( 2 ) (4»3) 


° may 

3ii order to find out p( 21 ) one/use the method suggested by Hill and Hunter 
(24) for single response case, using non informative prior distribution forO^ • 
Qn<»e p ) is obtained from equation 4 *3 a criterion similar to Box and 

Hill^''^’"^^^ may be arrived for multiresponse case with unknown variance 


covariance matrix. 
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SjBfTC HAIN NOPRNT 

C MAIN PROGRAM FOR I.IKEAR REGRESSION IN MUlTIRESPONSE SITUATION 
DlM£NSiON '63.C(2&)»3l 

COMMON"P»Y“(0O),X(5C»4 )/AREA/N,C<3,2§) ,eTA{2,25)/AREAl/S{50,50) 

1 

COMM0N/AReA2/SIGl(2i2),SIQ2(2,2),SI63'{2#2),SIG4(2»2) 

integer R,P,U,UU»UU1 

u ' EXPT BUN NO 

Q(1,U) OCKC CF A 

C(2fU) CCNC OF B 

C(3,U) CCNC OF 0 

NO OF RESPONSES 

ETA(l»Ur' rate of B 

ETA(2»Uj f^ATE CF 0 

N tOTAl i^o experimental runs 

P No OF PARAMETERS IN MQBELS 

Ra2 

PRINT OQS 

003 FORMAT(lHl) 

K«0 

READ 00i»IK 
C 

10 READ QOliN 
001 FORMATCfgy 

lF{N,eQ,63i PRINTC03 

K*K*1 

N2«2tN , 

C0NST*1,/F10AT(N) 

READ2',(Tc5<j;u>,^fl,3)i(|TA<I,U>,I?l#R)#U*l»N) 
prints ,aCU,UhU’^l#31i,<ETACI»y),I#l|RhU«l#Nj 
2 formatisF-t^s) ' 

PRINT' 0(jJ^;K' ’ 

004 rORMAT(lHO////49X,3iHLlNEAR REGRESSION WITH DATA SST,I3/49X, 34 tlH* 

1 ) ) ■ ■ . - 

C 

DO 80 M«l»4 
DO 20 I»1;N2 

DO 20 U»i;N2' 

sa»u)=o; ■ ■ , . 

IF(I,£Q',J) S(I,U)«1. 

20 CONTINUE" " ■ 

CALL'FORM(M,i) 

call. uinREginj 

SHliD'O. 

DO 40 U*i»N 
Y1*0> 

DO '30 U«1?P 
Y1»Y1*X(C}» 

3& cPNTiNuE" 

E10(U)»Y (Ul-Yl . 

SHl#l>sSl[(lii>4ElCtU)*E10<U> 

40 continue 

Sni»l)iOONST#SIClil) 

CALL F6RM(>t»§> 

CALL LIN'Rfc'GiN? 

SI (if2)=ur' ■ 

SI(2»2?'fO, 

DO 60 U®1»N 
Yl»0t 

DO 50 ^»l#P 
Yl»YltX<'U#U)#8{U> 

bQ continue" 

E20«Y(y>**Yl 



SU2,2)sSlC2»2)4e2C*E20 
60 CONTINUe""' 

SHi»2>»CONST»SI(l,2? 

SU2ti)»SUif2)' 

SU2*2i«5‘0{^ST*3l{2»2> 

PRINT 5g5Vr(S!( #J=l»2)f Is3l,2) 

PUNCH 5^5, USX(1,U)» Js 1,2), 1*1,2) 

555 F0RMAT<4^12,5) 

444 rORNAT?2§(2H«f ) ) 

C SI eSTjHATei)’ VAR-COVAR MATRIX OF OBS VECTOR IN EACH RUN 

call lNVSrsi,2) ■ ' 

GO TO' {6i;62,$3,64) ,M 

61 DO 66 Ii=l,2 
DO 66 01*1,2 

66 SIG3(Ii,ji>*SI(Il,Ul> 

GO TO 65" ■ ' 

62 DO 67 12=1,2 
DO 67 02*1,2 

67 SIG1U2,U2)»SI(I2,02) 

GO TO 65' 

63 DO 68 13*1,2 
DO 68 j3*i,2 

68 SIG2(U,43)iSHI3,U3) 

GO TO 65’ ' 

64 DO 69 14=1,2 
DO 69 J4*U2 

69 SIG4{I4',U4)sSI(I4,U4) 

65 PRINT444" ' ' 

00 70 Uil,N 
UUs2*U*l~ ■ 

UUl=UU*i 

s(uu,uu)»sia,i) 

. S(uui#uyj»sl|2ii) 
s<uu,uuiMShii2) 

S(UUi#UUlJf*^U2,2) 

70 CONTINUE 

0 ■ " 

CAUL F0RM<M,0) 

CALU lfNNEfe(N2) 

PRINT 007r i6tI),I*l,P,) 

C07 format <//,S0X,3x,4F12,4) 

PUNCH 555,«BU),IsJ,P) 

60 continue 

inK/3#3,£Q,K) PRINT 003 
CALL MATR1X(N,C) ’ 
inK,EQ,IKrC0"T0 900 
GO TO 10" 

960 STOP 
END 

SIBFTC FORM NOPRNT 

subroutinE'T&Rm<h,u 

C THIS SUB F0RMS'"THE OBS VECTOR YANDDESIGN MATRIX X FOR H.ppl H 
COMMON Pi V (50 J,X <50, 4 )/aREA/N,C(3 ,25) ,ETA'(2,25) 

integer P,U;UN 

psM*(6'.Mj/4' ' • 

IF (I,EQ,’0) GO TO 100 
GO TO il0#70),I 
10 DO 2ij U=l, N' 

Y(U)*ETA(1 ,U) 

X(U,1)*C(1,UJ 

X<U,2)s»'0(2i0) 

20 CONTINUE"' 

IF (M..E0,1) RETURN 
IF (M,8t,2j bQTO 40 
DO 30 0*1, N " 

x{u,3)*df3»y) 

3u CONTINUE 



RETURN 

40 DO SO UsliN 

XCU#3)i-»C(2»U) 

SO CONTI NtiE 

IF (H,lQ,3> RETURN 
DO 60 Uslf'N 

X(U,4)«C(3,y) 

60 CONTINUE' 

RETURN 

70 pap^(M*l>/2 

DO 8 0 U=ijtN 
X<U»1)*C12#U) 
au CONTINUE 

IF (M/2*2,NE,H) RETURN 
DO 90 

X(Ui2>s-CC3,U) 

90 CONTINUE-'* 

RETURN' 

ICO DO 110U»i»N 
UNpU+N ■*' 

Y<U)aETA(l»U) 

Y(UN)'«6TA{2iO) 
x(u»i)«cU,U) ■ 
x(UNiij'*o; 

X(U»2)a*C(2fU) 
ilO CONTINUE**’ '■ 

IF '{M;5T»2) GO TO 140 
DO 120 U»i»N 
UNaU*N ■” 

XCUNf2l»C(2,U> 

120 CONTINOE ’ ' 

IF (M.EQ.I) return 

00 ISO'D-i^N 

UN»U+N * ' 

x{u#3a»'Cf3,yi 

x(UN»3>'»-CcS;u) 

130 CONTINOE 
RETUf^N*" 

140 DO iSO U»1,N 

UN»U*N ’ " ■ 

X(UN,2)*0i 

X(U,3)aWC(2,U) 

X(UN|3)sCl^iU) 

15(3 CONTINUE-" 

IF <M,£Q,3) RETURN 
DO 160 U«iiN 
UNaU+N ' ■ 

X(U,4)8C(3#U) 

X(UN*4)a-0i3iU) 

16U CONTINUE 

return 

en d 

sibftc uinreg noprnt 

SUBROUTINE CINREGIN) 

C FOR UINEaR regression WHICH RETURNS THE VECTOR Of f EST PARA 

C METERS 

COMMON PiY(50) »X(50i4 I/AREAl/'St 50|50) ,B(4 ) 

INTEGER p‘ * ' 

DI ME NS lO N XT U >5 0 ) f X TS (4 #5 0 ) » X TS X( 4, 4) #x TS ti 4) 

DO 10 lalfp” 

DO 10 U*i»N 
xT(iij?axu,n 

10 continOs ' 

C S ' iStlRSi OF VAR COVAR MATRIX Of OfS VICTOR T 
DO 20 I^llP 



PO 20 K*1#N 

XT5(l»»^)»XTS(IfW)*XTCI»K)*S<K,j) 

20 continue' 

DO 30 lEifp 
pO 30 ijsi t P 
XTSXdiWfsO, 

DO 30' KEi,N 

XTSX(I» J?«XTSXd#U)+XTS(I,K?*X(K» j) 

30 CONTiNUis " . 

CALL INVy{XTSX,P) 

DO 4 O' I'a'iiP 
XTSY<I)aC, 

DO 40 J=1.N 

XTSY(UsXTSY( I)4.XTSU»J)*Y<U) 

40 continue' 

DO 50 lsl»P 
BCI)«0', " 

DO 50 U*i|P 

8Cn=BUHl)<TSX<I,^I#XTSt(U) 

SO CONTINQE "" 

RETURN' ■ 

END 

SIBFTC INVS NOPRNT 

SUBROUTINE’ iNVS{A,M) 

C SUB FOS 'lRvS"0f 'MAt 'A 

Dimension A(4, 4)' 

DO 60 KsiiH"" ' 

A(K#K)i»»i;/A(K#K) 

DO 20 ■■ 

IF{ I-K)'1'!J#20#10 
10 AC IiK)»*A(tpK5tA(KiX) 

20 CONTlNilB 

DO 40 iPliM 
DO 40 u4ii'H 

IF(Cl*|(?tj>KI) 30»40#30 
30 AC|,U)fAU#J?»A(!|K)«A{K»U) 

40 coNTlfiUe ' ' ■ ■ • ' 

DO 60 JEIiM 
IF(U-K558#60i90 
50 A(KfU)E-A(IC, J)#ACK»K) 

6U CONTINUE' ■ "" ' ' 

DO 70 1=1, M 
DO 70 

A C I I J ? a A { 1 1 U ) 

70 CONTINUE 

RETURN 

END 

SIBFTC MATRIX NOPRNT 

SUBROUTINE HATRIXCNiX) 

DIMENSl0f!!''X(3,N),Xl4C25,3)#Xl2C25,3),XTilCS»25},XT12c3,25l, 

1 AMXll(3,a) -, AMX12C3,3) # AMX21CS,3) ; AHX22C3#3I ’ MAlc3i 3) 

2 , X21(25|3l , X22C25#3),XT2l{3,2S)#XT22Cl4SJ#iHXllC3|3r ' 

3 , 8MXi2c3",S) , 8MX2lc3»3);8MXg2(3,3),HA2(3,3I,X31c25i2liX32c25,2) 

4 , XT31<2,25),XT32C2,25),CMXllC2#2)tC«Xl2{2,2J,iSMX2iC2,2j# 

5 CMX22(2#2),MA3(2,2j,X4i(25,4S,X42i2f|4),Xt4iC4,25j;xT42C4,25) 
DIMENSION dMXai(4,4),DHXi2<4f4) ,DMX21 {4, 4 1 , DMX22(4;4) ifHAA (4,4) , Ai,l 

X (3,3),BLMf»3),‘CLi(3,3),Da{3,3),AU2{3,3),Bl2(3,3)fCL2(3;3);DL2<3 

2 ,3)#AL3(2i2),BU3'<2f2),CU3(2*2),DL3'(2,2),AL4(4i4),8L4’(4;4),CL4C4,4 

3 ) ,DL4 (4^4 ),Bi(4,l) ' ' ' ' ' " ‘ 

C0MM0N/AREA2/SIB1(2,2),SIG2(2,2),SIG3(2,2),SIS4C2,2) 

REAL MAi#MA2,HA3,MA4' ' ' ' " ’ 

Mi=3 

M2»3 

M3«2 

C MAIN LINE PROGRAM 

WRl'Tft '6f !? )( CXCI >U),Ual.N) »i=1.3) 



97 


333 

999 


F0RHAT(15F6»2J 

CAll 

TRF6S£(XllfN»Hl|XTll) 

TRF0sE(Xi2iN,Kl,XTl2) 

MAMUi.CXtlliXli#N#Hl,AMXll) 

MAM0i.«XTii,Xl2,N,Mi,AHXi23 
MAMUl(XTi2,Xll,N,Ml|AMX21) 

MA«ULCXta£,Xa.2,X-,Mi, AMX22) 

TRP6at(X21#N,H2,XT2J.) ' ’ 

TRPpS|cX22,N,K2,XT22) 

MAMyi.<Xt2l,X21»N»M2,8MXll) 

MAMyUXt2i,X22,N,M2,8MXl2) 

MAHyUXT22#X2a#NiM2reMX2i) 

MAM0LtXT22,X22»N|M2|6MX22) 

Mr uKMt 0# X# Nf X3 1# X32j M3 ) 

TSP0SE(X31,k,m3,XT31) 

Tftp0s|(X3£,N,HS,xT3a} 

MAMUC(Xt3i,X3ifN>N3,CMXll) 

MAM0ltXt3l#X32,N,HS,CMXi2J 
MAMyi;cXT32,X31,N,M3»CMX2i) 

MAM0llXt32,X32,NfH3iCMX22) 

^Ay?CftXii( CMXi2f CHX2iifCMX22fS}G3fH3fMA3f Ai,3 f 6|^3 f CL3 « Ci,»3 ) 
MF0^M(4»XfN,XAi,X42,M4) ■ - 

TRP0Sfe(X4l#N»M4,XT4i3~ ' 

TRP0s 6(X42#N»H4,XT42 j 

CALI, MaHUC C XT? i i X4i; N f M4 , DHXil ) 

CALL MAM0L{XT4i,X42#N,M4,fiHXia3 
MAHDaxt^2#X4i,N,M4,CiMX2iJ 
MAM0l-(XT42fX42#NiH4,DHX22) 

MAT|UflXi|»DMX12»DHX2i,DMX22,Sj64,M4rMA4,AL4,BL4,CL4#Dl.4) 

MAtiNV'(MAl,3,0i,O.DiTER ) ■ ... - • - 

MAT INV<NA2,3| 0i,(3,|3efER) 

MAtrNViMA3,2,S3„c#D£TER> 

_ HAtTNVtHA4,4,Si,9,DETEft) 

WRlT|(6,3j'3r (■(MAlU,y);4=lf33,I»l,3) 

WRlTEUrSSs) UMA2(l,Jh0ni3M*l|3) 
wRiTE(6,3i3) 

WRITE (6# 333) ((MA4.Uf;jJf j4i#4hl«l,4) 

PUNCH 333,UMAi<r,4'J»j«i;3)aPi#3) ’ 

PUNCH 333# crMA2(T»U?#0*l,3),’l = i,3) 

PUNCH 333,((HA3a#;5)i j*l#2),tPi#2) 

PUNCH 333#C.{HA4(T»U)»U=1,4)#I*1,4J 
F0RMAT<6E12;5)’ ■ * ' ' • ■ ■ ■ ■ 

WRIT6{6"#9$9) ' ' 

FORMAt(58'<2H##) ) 

RETURN' 

end 


CALL 

CALL 

CALL 

Call 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

call 

CALL 

call 

CALL 

call 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

call 

CALL 

CALt, 

CALL 


call 

CALL 

call 

CAU 

call 

CALI. 

CALL 


C 

ilBFTC MFORH NOPRNT 

SUBROUTINE’ !ilPORM(IMODELf X#N#2l#22,N) 

dimension, X(3,N),|ia.<N,M),22tN,H) ' 
IFUMODIL.nI’.I) Gp TO 'lOO ' ' ■ 

Do 3,0 U#i#N 
Z1 (U#.3,)**^tlid) 

21(U#2)s0f'‘ 

Z 1 g#3 )»Q, 

Z2(U#1)'«X<1#U) 

Z2(U#2S»-X(2;0) 

10 Z2{U#3)iix’(|#u| 

iOO IF( lHpOI|;*N6i2) CO TO 101 
DO 3L| Uii»N " 
zi{j»i)»-xa#j) 

ZlU»2j»xU'tJ) 

Z1U#3)»0, 



Z2( 

22{j»2Js-:K(2j j) 

11 Z2(J#.3?«^)(Uf0j 

101 IF UM50&1;ne;3) go to 102 
DO 12 Js'i^N 
21<JtlJs-X(l»J) 

2i<J,2?sO, 

Z2<>)#l?sX(l,0) 

Z2CJ,2j«-J{<2,a) 

12 CONTINUE 

102 IF <IM0U£L,NE,4) GO TO 103 
DO 13 jsi,N 

2i(j#i)s-xa»u) 

Z1(Ji2}5X(2;uJ 

zi(j#3j*{j;' 

Zl(j#4>s0, 

Z2{U|l)sX(l#J) 

Z2{U»25»-X(2JU) 

ZZU^SM’-'XUii) 

13 Z2(U#4JEX’C3;U) 

103 return • ’ ' 

END 


c 

IIBFTC TRPOSE NOPRNT 

SUBROUtlNfc tRP0S6(lNMAT» iNjIHfOUTMAT) 
dimension" !NfiAT'UN/lH^,OUTMAtCIH»IN> 
REAL INMAf ' ' ’ 

DO 10 IbIHN 
DO 10 Qi'iftM 

10 OUTMAT<Uny»INHAT{I»U) 

REtURS 

ENO 

C 

SilBFTC MAMUL NOPBNT 

SUBROUTINE' flAHUL{AfB,lN» IM,C) 
DIMENSION A'(lii,INj,6(IN,lN>#C{IM,IM) 

DO 10 iRliIA 

DO 10 U«i»'lH 

c(uu)=o; "■ 

DO 10 'K*l# IN 

10 C(Ifvl)a.dU'fU)+AUf«)*P<K»U) 

RETURN ■ ■ ■ ' • • ' • * ■ 

END 


C 

SIBFTC MAT NOPRNT 

SUBROUTINE" MAT(ArB»C,DtSISMA,MiMAiAL#8L#0L#DL) 

DlMENB|0N''A(M;M)»0lM,M)»d<M,M)»D(H;Nj,AL(MiH>»8L(H,MI,CL(M,,M),DL{M 
1,M) iSISMA{2;2),MA(M»M) ■ ■ ■ • • 

REAL MA ' 

DO 10 I»1»K 
DO 10 U»i»M 

AL (I W)*SlGMA(lil) *Aa tJ) 

BL(I,Ul»S'lgHA(i#2)*f n#U) 

CL{I»Uj-^S’lffMAUilHd(I#U>’ 

DLUiUi«|XSMA{2»2M&il»^> 

10 MA(iiU ) * AL (IfUl+BLdfUI^CLtifU J♦Dl» ( I # U ) 

RETURN" ' ■ 

End 


SIBFTC MATINV NQPRNT . „ , 

‘ SUBROUTINE' HATINV{A#N,e#M#pT|RMI 

DIHENIION A(N|N)»l(N#MhJPiVpT(40)f INDEX(40#2) 
EQUlVAC|N'UE‘a'ROH,^RON), (AHAX,|*fHAtl 

10 dETERMiI#?' 

15 DO 20 4»t»N 
20 IPIVCT(wj = U 

C SEARCH fUH PIVCT El.&^E(^T 
30 DO 5E50 



40 

45 

50 

60 

70 

80 

85 

90 

95 

100 

105 

110 

130 

140 

150 

160 

170 

200 

205 

210 

220 

230 

250 

260 

270 

310 

32o 

330 

34o 

3»0 

355 

360 

3^0 

38 0 

39 0 
400 
420 
430 
450 
455 
460 
500 
550 

600 
610 
620 
63 0 
650 
660 
670 
700 
705 
710 


11 

12 

991 

m 


SENTRY 


AMAX»0,0 
DO 105 jal#N 

IF(iP|VOf?^)-»l) 60,105 ,60 
DO 101 K«i,’N 

IF(IPIV0T(K)-1) 80, 100, 740 
fFJAHAX-ABsUUiK))) 85,l0O, IQO 
IROW«J 
IC01UH=K 

AMAX=ASS<A(g,K)) 

CONTINUE 
CONT INUE 

IPIVOTUupLUH )sIPIVOT{ICOlUH )♦! 
interchange RONS TO PUt' 'PlVQT ElEHENT 
IFC iRON-iCbUUM } 140, 260,140 
DETERMs-mTERH ' ■ 

DO 200 Nsi^N"' 

SWAP«A(IR$N,U) 

A(IROW,'U’A(ICOUUH , 1 ) 

AUCOUUH ,irp§NAP 
IF(M} 260, '2601210 
DO 250 

SWAP»B( lR6N,t) 

B<lRdW,i,?»BMCOlUN ,1.) 
B<ICOUUM',Lj»§NAP ■ 

INDEXn,l)'»ii^ON 

INDEXtIiaj«''ieoLUH 

DividE' PHdr"Rdw iY pivot elehent 

PIVOT »AnCOlOH ildOLUM ) 

OETRN»bitpM'#p'iV0T 
A{ICOl‘^M?|Cbl.0M)Ai,C 
DO 350 

AUCOtUH»'C)*A<ICCUUM,l.^/PIVOT 

IFCM) IBH, 360,'36‘o‘ 

DO STO'C^lift" 

8(l60tUH.,'C)!=BUC0L.UH,l.)/ PIVOT 

REDUCE NOK 'PIVOT' RONS' 

DO SSQ 

IFCUI-IOOCIJ^) 40 0, 550 ,400 
TeACU,lcOCUrtJ 
Aaif ICO'CUH|aO,0 
DO 450’i:*i,N' 

A<Ll,U)sA^Ul»U)•A(ICOLUM,U)4T 
IF(M) 55o;''550, 46r ' ' ' ' ' ' 

DO 500 ’C*1,M ■ ■ 

B(Ui^U)»8IUi,l)-S(lCOl.UM,l)*T 

CONTINUE ' • 

interchange COUUMS 

DO 710 ■ I»r,'N 

i«N*i-r' " ' 

IFaNDexa,l>-iNDEX(U,2)) 630, 710, 630 

jR0W=?lNDEXct;,2) 

DO 705 Kii;K ' ' 

SWAP*A(K,'<3rON) 

A(K,aROWI«A(K,UCOUUN) 

ACKfUCOLOhliSWAP ' " 
continue" 

CONTINUE 
DO 11 k*liN 

IFdPIVOtIK) ,NE,1) SO TO 12 
CONTINOS" ' ' 

RETURN” " 

PRINT 9fl 

rORHAT|/IOX»MA|RlX I.l SINGULAR#/) 

return ~ * 


ON DIASONAI 



« o 


SWaTfOr M6F287 
IBjOB 

IBfTC MAIN 

MAIN IINB PROGRAM FOR COMPIKX AISORITHM OF BOX 
dimension X(6;4>,R{6,3),F{<S)rG(4),H{4hXC<3l 
COMMON TH|tAM»4},P(4),0SI(4i4>,DETi»DET2,DET3»DET4 

1 “ .M3,3>,B(2,2),C<3,3l,D{2.2),|'(a,2)rA 

2 AH ( 2 f 2 ) 

INTEGER GAMMA 
Nl=3 

N2s3 
N35 2 
N4»4 
!'^5=4 
N6b4 
N7=2 

NM0DEU»4 


NI«5 
NO* 6 


READ(5 4)N»M,K,ITMAX,IC»IPRINT 

1 FORHATIBiS)' ■ ■ . • ■ 

REA0I5' ;2)' AUPHa^BETA, GAMMA 

2 FORMATISiElO,'4,l5)'' ‘ 

DektA*e;dH " ' 

4 F0RMAT(8Ei0.4) 

DO l01"it«2.K" 

REaD( 5 4r(R(U,g4,^j«l,N) 

3 F0RMAT{3F8.4) ' ■ 

101 CONTINOr' ■" 

DO 3^9 IDIF*l,9 

mf 

104 FORMAtUK,«DAYA’ SIT, #,I2) 

READ?$'»10O)f <A(l4)#'j*l,Nl)#I»l,Nl) 

REaD(5 400) a@U4)4»l,N7),I»lfN7) 

READ(5 »|0O)(<C(J4)4*1iN2)/I!*'3.,n2) 

READ <5 ■ 400 )UD {14)4447) *1447) 

READ (5 *|00) {(S{J»4»j*3-43)*Isi*N3) 

REaDIS 400) ({AF'<I4)4«i*N7),l»l,N7) 
rEAD< 5400) t {AG'(It^)4*i44),I*l,N4> 
READ<5400){<AHUf J)*a=l*N7)*I«l,N7) 

REaD( 5 *iqo)(P(l)»iniNMODE4‘ 

READ'tSjloE) (4 HETA(i4)|J«1,4)#I«1#4) 

100 F0RMAT(6Ei24J 

102 F0RNAt(4|i24) 

•«y . 

DO 299 ITER *1*2 
READIS ,4 (X{i,W)*%l*l»N) 

WRITE(6,3:03)'(X(1»4'*U4*N) 

103 F0RMAT(2X*4lNlTlAL X ARE *’ ,3(2x*|J24)) 

NRITE'ce ,10) 

10 formats * ' 18X,24MC0MPLEX PROCEDURE OF lOX) 
WRITE{S’*18) 

18 formats ■■ 2X40HPARAMETERS) 

WRlTE(6’,ll)N44,ITMAX,lC»Ai.PHA*BiTA,eAMMA,DElTA 

11 formats '"SX^HN »",I2#3X,4HM' ? *!2,?X*4HK ' • ' * 124 

1 14,2X4H1C • ,U,/42X#8HALPH^ * ,FS»2#SX,7HiEYA ^ 
aSHOAMMA '»”*I2i3X,8HDEUTA * 48*8) 

lF(lPRINt)4d»50»40 
40 WRITES'^ 42) 

12 F0RMAtS/42X*14HRAND0M NDMBiRS) 

DO ,800' 42*IS ' ' . 

WRlTe(6 »i3)(U*IiRCUiU#I»l*N) 

13 FCRHAT(/4S2X,2hF(42»lH* .12.4H) * ,F8,4»2X)) 

800 CONTINUE 

0 


F(2*24A§<4*4) 


,8RITHAX * , 
tFio4,3x* ■ 



o Cl Cl rs o o o’O ci'Ci o-o ci 


50 CAUk C0NSX<N,M#K# nHAXfAl,PHA|BETA»fiAMHAfDfclTA^XtR*f 
1N0»G#H#XC,JPRINT)' ■ ' 

IF( 2C,2C»S0 

14 ’^FORHAT{''^^'X#^ChnNAl. VALUE OP THE FUNCTION ■ ,E2Q,8) 

15 ^FORMAT? ^ ' SXilAHFIlvAL X VALUES) 

DO 300 J»i»N 
NRnE(6'#l^> U»X(lfcV2,J) 

16 F0RMAT</» 2Xi2hX (» I2|4H) « »E2(5,8) 

300 CONTiNOt 

CAUL PUHRQ<X,16V2#1?IF) 

GO TO 99? 

^ .. . 

17 ‘'^FORMAn/^J^ExIsSHTHE NUMBER OF ITERATIONS HAS EXCEEDED iH^lOX# 

1 10HPROORAH TERMINATED) 

999 CONTINUE' 

WRITE(6|106) IT 

106 FORMAT(2X,flTERATI0N NO «*#H) 

299 CONTINUE ’ 

399 CONTINUE 
STOP 
END 

$Ii;iFTC^C0N|X^^^^ CpNSX(N»M#Ki ITMAXi ALPHA t8ETAfGAMMA»DELTA?X»RrF? IT# 

^COORDINAIIi’III&IAI^PURPOSE SUBROUTINES 
AROUEMENT LIST 

lEVl ’*'»'!w1x^0F^P0INT with MINIMUM FUNCTION VALUE 
IEV2 « i>?EX OF POINT Wlt^ maximum FUNCtlfaN VALUE 

KODE ^^CONTROL^MEY USED TP DETERMINE If IMPLUCI.T CONSTRAINTS 

ARE PROVlUeU’ 

K1 ■ Up LOOP LIMIT 

all other PREVIOUSLY DEFINED IN MAIN LINE 

DIMENSION Xl6i4)»Rt6»3)#F (6) |G(4) »HC A) #XC{3) iXBESTlS) 

integer gamma 

u 

n=i 

K0DE=0 

IF<M*N)2Q»204Q 
10 KODEsi' 

20 CONTINUE 
’ DO 40 

DO 30 J'sl^N 

3 0 X( II »U>A0*0 

4 0 CONTINUE 

CALCUU^Tt COMPLfeX POINTS ANO CHECK ASAINST CONSTRAINT? 

DO 65 IIPZ.K 

DO 5 0 U'lfN 

1* 9 f f 

GALL CONST(N»MiK#X#GrHfn . 

XUl»U>A§IU)4 R<ll#U)«(HIJ)*GtU>) 

50 CONTINUE 

call Check <n ,m ,k ,x ,s , i ,K0DE» xc ,D gLTA»«ff ^ 

IF C 1 1- 2) 51 is 1, 55 


0 

c 

c 



o Oci C'-OO 


18 FORMAT(//»2X»2CHCOCRDINATES OF INITIAL COMPLEX! 

lOSil " ’ 

WRITE($ ,19!(IC»>i»XCI0,*J)#Usl»N) 

19 F0RMAt{/»2(2X,2hX<#l2,iH## 12,4H) » ilPEi3*$)! 

$5 IFCIPR|NT)E>6,65,66' ' . - 

H WRITEli& II?) J-1»N) 

6S CONTINUE ' ..... 

Kl^K' 

DO 70 I=1»K 

CALL FUNC{N»M»K»X»F#n 
7D continue 
KOUNTsl' 

IA=0 

FIND POINT NITM LON&ST FUNCTION VALUE 

IFUPRINT) 72»8y,72 
72 WRITE(6 ,'2lj’ 

21 F0RHAIU#2X,22HVAIUES OF THE FUNCTION) 
WRITE<6',22)<u,F{U)#U»1#K) 

22 F0RMAT(/#2<2X,2HF<,12,4H) » ,lPEi3.6)) 

80 IEVl«l ' 

DO 100 ICM»2,K 

IF<F(ieV|)'FUCM)) 100»100 j90 

90 lEVlsICM ■ 

100 CONtlNUE 

FIND POINT WITH HIOHEST FUNCTION VALUE 
IEV2«1 

DO m ICM«2»K 

IF(F{IEV2)-F(I‘CM)) 11C»110i12 0 

110 IEV2»ICM' ' ■ 

i$o continue 

DO iiigQ«!i,N 

111 XBEST(JuFEX<IEV2,UU) 

C CHECK CO NVER’OInCE" criteria 

IF<F<iev2>-(FnEVl)t8ETA?) 140,130#130 

130 KOUNT»i 

60 TO ISO 
140 kount«kOunt^i 

IF(KOUNT"CAHMA)150»240i240 
C REPLACE>SiNT WITH LOWEST FUNCTION VALUE 

150 CALL CtNfR(N;M,KiIEVl,I,XC,X»Kl) 

DO 160 43®i#N 

160 X(lEVl,ODI«Ci,0+ALPHA)tCXC(UU))» ALPHA* (XdEVlfUU)) 

idEVl ■ ' ■ 

CALL CHECK IN# M#K# X.6^H» If KQDEiXC, DELTA# Kl) 

CALL FUNC(N#K#K#X#F#I) ■ ' ' ' ' ' 

C replace N8W POINT IF' IT REPEATS AS LOWEST FUNCTION VALUE 

ICHEKs'O " 

IKNT*0 
170 IEV2al 

ICHEKs«lGHEK+l 
DO 190"lCB52iK 

inFaEV2')-R<lCM))190#190#lSO 
180 ieV2^iICM‘ 

19 0 CONTINUE 

IFU eV2iJfcVl) 22 0, 20 0# 22 0 
200 DO '2lO UU'«liN 

X(IEVl»U^)‘»CXUEVl#UU)*XC(UU))/2,0 

210 continue 

|ii|f Vl" 

Utl !HlCK(N,M,K,X,l,H,I,KOPIiXC,DflTA,^l) 
call fu;\c(N»m#KiX#f. d 
IF(lCHtK-25) 170,202,202 
202 IKNT«IKNT*l' 



WRlTE{6i2U3? IKNT 
203 F0RMAT(2X#i|KNT* ##I3) 

peilRN 

00 410 jC *1»N 

X(IEVl|0i:> atXUEVliOtJ +XBEST(^L))/2,a 
410 CONTiNUfc' 

GO TO 170 
220 CONTINUE 

IF<IPRINT) 230#228,23Ci 
230 PRINT UgSin 
23 FORMAKU) 

PRINT ij24’ 

C24 FCRMAT{/72X,3ghC0CR0lNATES OF CORRECTED POINT) 
print 01?»'C'lSVl#UC»X(IEVl,UC),JCsa,N) 

PRINT 021’ ■ ' 

PRINT 022i(I»Fa),in#K) 

PRINT 02^ " ■ 

025 F0RMAT(/#2X,27HC00RDINATES OF THE CENTROID) 

PRINT 02ei{UC,XC(UC)iWCn,N')" 

026 F0RMAT(/#3{55?,2fiX(f I2;6H#C) * ,E14,6#4X)) 

228 IT*IT*i ■ ■ ■ 

IF(n-lTHAX)80,80|24D 

240 RETURN " 

END 

sibftc check 

SUBROUTINE CHECK(NfH,K?X,6»H,ItK0DE|XCiDElTA,Ki) 
C ARQUEHl'NfolSt • 

C ALL ARGUEMENTS DEFINED IN HAIN LINE AND CONX 
DlMENSlO?^' X(§j4hB'(4),H(4),XC(3) 

10 Kt»0 

CALL CONST(N,M,K#XrG»H#I) 

C CHECK ACaInST' EXPLICIT CONSTRAINTS 

DO 50 UiiVN'"'' 

IF{XaiJ)*QU)) 2C»2Q*30 
20 XdisIJisGiUHOELTA 
GO TO 50" ■ 

30 IF(H(U)-XU»U) )40#40,50 
40 X(IivJ)9H'<Oj-0EltA 
50 CONTINUE" ' 

IF(KODS)110»il0#4O 

C CHECK AGAINST' IKPLICIT CONSTRIANTS 

60 NN«5N*1 

DO 100 U«NN»M 
CALL CONSt(N#MrK#XrGfHiI) 
iF<x<i,g)SQ(a))ee#7e,7o 
70 lF(H«U)»)^{|»^h 8C#100»100 

60 lEVl^t 

KT*1 

CALL CENTR(NiM#K#IEVl#IfXC,X#Kl) 

DO 90 UU'ijN'’ ' ' 

x(iiUu56(xa»uu)+xc(ju))/2,o 
9 0 CO NT IN UE 

10 0 continue 

IF«KT)iiO»iaO»10 

110 return 

END 

sibftc centr 

subroutine CENTR{N#H,K,IEVltI»XC,X,Kl) 
DIM6NSI0 n"x]Si4>iXCU) 

DO 20 Ua#N ■■ ‘ 

XCCU)*||0 
DO 10 IL»1»K| 

20 XCtU)=(XCI J)-X<IEV1»U) )/(RK*-l*OJ 
RETURN 
END 



*|oFTC FUKC 

SUBROUTINE FUNC(N,H,K,X,r,n 

DlMfcNSiCN~X(6#4),F{6) 

COMMON THETA(4#4JiP( 4),0SI(4>4) ,DETlrDET2,DET3» 0ET4# 

1 u,o O. A<3f3>»B(2,2),C(3,3),D(2,2),E{2,2),Ar<2,2), Q(4,4), 

2 H (2#2) ’ • . - . . ^ • 

xi«xafi) 

X2*X(Ii2j 

X3«X<I,3J 

xii»xi*xi 

X22«X2#X2 

X33aX3#X3 

X12»X1*X2 

X23»X2*x3 

X13«(X1«X3 

C ELEHENtS OF VARIANCE COVARIANCE MATRIX OF MODEL IFOR (N4i)RUN 
Asui»§aii)+A(i,i)#xu' ■ . ... 

ASU2»B(1 |2)*A<a»l)#XiitA<l#2)#X12*A{ii3)*Xl3 
ASI2X>eASliE ■ " • • • • 

AS|22»E(2,'a)4A(l,l)#XH*AI2,2)#X22iA(3,3?tX33*2,#ACl#2)tX12 

l«2,#A(2f3)#X23+2,#A(l,S)#Xl3 

C elements or variance covariance MATRIXOF MODEL 2 FOR (N+ilRMN 

BSUl«0(lfl)+0(i,l)»XU*C«2,2>#X22«2;#cU»2)*Xi2~ 

8SU2»D^»2j*g(i,i)#xil*IC(2#2)*d(2|3^)4X22i(2^#CCl|2)*C(l#3) )»X12 

BSiais'eSiiE ' " ’ ■ ' ' ' 

BSl2gAD'(2;2)tC(l,l)#Xil+(C(2,2)*2*»C<2,3)^Ct3,3} )*X22»2,»<Ca»2H 

1 C{i,3))#xia 

C ELEMlNfS' OF VARIANCE COVARIANCE MATRIX OF MODEL3FOR(N*i)RUN 
CSlURlUVl)#XliAAFU,i) ' ’ 

CSU?»AFfi#2)*Ea4J*Xll*E(l»2)#Xl2 

csiaiAcsila ' ; ' ■ ' " ' 

Csla2AAF’(i»2)*E(l#i)*Xli»2i#E(l»2)<iX12*E(2»2)#X22 
C ELEMENTS'OF VaRIANCE covariance matrix of mODEI,4FOR(N*1I run 

DSUl*rtU#l)+Sci,l)*Xli*2i*G(if2)*xl2+S(2t2’)#)i(22 ; 

DS U2»H( i» 2) -S (1 4 )*xai* (2 >61 U 2) ♦O (1 ,3 )) *Xi2»( St 2# 3) *Q (2 #2 ) ) »%22 

l«G(i,4j#X3.3tOt2f4)«X23 

DS 121*0$ ‘I'i2' 

DSl22»!H<2»2)*G(l,l)<»Xlit (G{2,2)*2«tS(2,3)*e(3,3) )#X22*Gt4, 4? ♦X33»2 

l,*{Gt2;4j*5t3f4 4#x23*2,*G{l,4)<»X13 -2,0#(GU«2jAO(i#3n»Xi2 

C DETERMINANf OF’ THE' MATRIX ' " ' 

DeTl*A$lJiiASl22»ASU2#Asil2 
DET2sBS'lil#BSia2»BSU24B$Il2 
DET3*CsUi»0sr2a-C$U2#CsU2 
D ET4= DB lii* 6$ l’22»DS 112403112 

C ELEMENTS’OFINVERSES of variance covariance matrix of .MODEL! 

OASlt«ASl22/uETi 

OAS12swASU2;'D£T1 

OAS2l»OAS’a;2' 

OAS22sASUiXDETl 

C eLEMENTS'e? INVERSES OF VARIANCE COVARIANCE MATRIX OF M0DEL2 
0BS11«ISI22/DET2 ' " , , ~ - 

08 SI 2* *8 'Si 12/0 Et2 
OBS2laQBSii2' ■ ■ 

08S22*SSfii/BET2 

C ELEMENtS'Or INVERSES OF VARIANCE COVARIANCE MATRIX OF MODELS 
□ CSll»CSl22/iiET3 
OCS12««CSil2/DEt3 

00321*00 si 2' 

OCS22»CSlii/0BT3 

C ELEMENTS'OF INVERSES OF VARIANCE COVARIANCE MATRIX OF MODEL# 

ODSll*&SI22/OEt4 ■ ■ •• - 

ODSi2«-0'Sii2/ilf4 

ODS21«Of|P' 

0DS22*bSUl/DET4 
osi\i»i^“U As li 
OSI(li2)*:PrtSl2 
0Si<l*3) spAS 21 



.•:SI {l,4)aOAS22 
OSI (2,iJ«0bSli 
OSn2»2)iOBSi2 
OSI {2i3jab6S2l 
0SI(2#4j«b6S22 
osu3ii j*^bosii 
OSI(3#2ja0cSi2 
OSI (3,3 HcCs 21 
OSU3,4)»Oi5S22 
osn4j.iMbi3sii 
0Sn4#2)abl!Si2 
OSI {4,3 j*bus2l 
OSI {4#4)=bbS22 
C EXPRESS lO'N 'for Ya-Y£ 

Y1.2a-(TH|TAri#l)«7N|TAI2,l)) tXl-THETA (2 ,2 )*X2 

Y21«(THETSa»l>-THEtA{2,l))»Xi-(THETAa^2)-TH6TAC2,2)-THETA<2,3H<* 

1 X2*THeTA{lf 3)*X3 ' ' - ■ - 

C eXPRESSiOR for Y1-Y3 

Y13a-(tHfc;TA(l,l}-TN6TA(3,lI?#Xi 

Y3i«ITftEtA(i#l)*TH£tA(3,l) ji»XU{THETA{l*2)*THETA(3,2j )*X2*THETA (1 # 

13)#X3 ■ ' ' 

C EXPRESSION FOR Y1-Y4 

Y148«CTHEtAUil)*TliETA(4,l))#Xi-THETA(4,2)fX2 

Y41s(Tft£TA(i;i)*THEtA(4,l))»Xi-(THETAUt2)-THETA{4»2)-THETA(4,3)) 

l*X2*iTHEtAU»3)»tHEtA(4,4j j#X$ ’ ' ' 

C EXPRESS! blY'PoR Y2 -y3 ' ‘ 

Y23e-ItHETA<2»i)ATH|TA(3»l))*XitTHETA{2,2)fX2 

Y32*|TftETA(2a)»THEtA(3,l))«Xi«(THStAia»2)-tH£TA(3#2)*THETA{2,3))4 

i X2 - ■ ■ ■ • • ■ ■ • 

C EXPRESSION FOR Y2-Y4 

Y24*-«tHlTA?i2»l)-TH|TA(4,l))#Xl*i[TH|TA{2,2)-THETA(4i2))#xa 

Y42»{TR£|S(2a>-THETA<4#i)H^«i»(tHETA<2»2)«TRETAi4»2)-Tl^lEtA(4,3)-t*t 

i H6fMa»3n#X2«<THEtA<4,4j4X3) ' 

C eXPRlSSlONFOR Y3*Y4 ' 

Y34ii«l[fHElA(3Ji')-Tk|TA<4,l) J#X1«THETA{4|2HX2 

Y43#(Tfi6TAU#iWTRETA(4;a)}#Xi*cTHETAUr2JAtRETA(4»2>»TWeTA(4,3)># 

lX2-fHetA?4;4>#X3'" ' 

C EVAlUAtI3N”'0F' K¥ 

C FOR Pi ANO PS' ' ... 

TRAi2*ASll3.t0BSai+A8li2*0BS12tSSUl#0ASU*8Sll2#0AS12«2j, 
TRA2i5ASti2*o6S12+AS'l22#08S22*asU240AS12t8SY224o4S22**2t , , 

AMllaOASii+OBSli' 

AH12sCASi2+08Si2 

AM21*AH12 

AM22sOAS22^0BS22 

RES3.2»AMii»Yi2#Yl2*2,«AM12*Y12tY21tAM22#Y2i#Y2l 

P12!*P(a.)4P(2)#{tRAi2»TRA21+RE$i2) 

C FOR Pi Ai^0'P3' 

TRAl3s»ASlii*OCSai^ASU2»OCS12*CSU1^0ASll+CSU2«OAS3.2-2^ 

TRA31sASYi2«bCSl2*AsU2»0CS22*0sfl2*0A$i2+0sl22#OAS22-2, 

BMll«OASiitObSu ^ ' ' ” “ ' ’ 

BM12?OASl2*bC$12 

BH22«dAS2^AoOS22 

RESl388Mii»Yi3*Y13*2,#8M12^Yl3*Y3l*BM224Y31#Y3i 
P13»P(l)iP<3)4{TRAl34TRASl4RESl3 J ’ ' 

c FOR Pa . "and PA' 

TRAa.48ASmAOCSia4ASIl240DS12*DSIll#OASH4DSI12tOAS12'»2i 
TR A4 1* AS’ii 2* Ofi si 2* AS 12 2*» OD S2 24 pS II 2# OA Si 2* 12 2# OA $2 2* 2, 

CMil*OASl 1400511 
CM12«0ASi2400S12 
CMElaOHlS:"' 

CM22sOAS2240DS22 

RESl4«CNii4Vi44Yi442,*CMl2#Yl4fY4l4g«||l^^ii|..<,i 
Pl4»P( i)4P(4,«(tRAl4*TRA4i4R£fi4V ,, 

Ti ij.<»ocsi2+CBSlltCSIll+CSIl2tOBS12«2i 
^i43?.8BS Ua*0CSl2+^ Y<»OCS22-(-CSla.2«CBSl24Cs'|2S*0BS£S*-2* 

' iii, 11 ...... 



D^'^iisObsii+ocsia. 

Dl^l2*08Si2+C}CSl2 

DM21S0M12"' 

DH22»0BSa2+pCS22 

RfcS23aDMii*V2ifY22#2,«DH3.2<»Y23#Y32*DM22<»Y32tY32 
P23« P(2)*I^(3)»{TRA23*TRA32+RES23) 

C FOR P2 AN0’P4' ' ' 

TRA24a8Slil#0DSU+BSI12#ODS12+DSIli<»OBSU*DSU2«O8S12<-2!» 
TRA42»BSli2*O0Si2*BSl22»ODS22+PSU2*OSSi2*0$l22*d8S22-2, 
EHll«OfiSiltdOSll ' ' ■■ ■ • ■ 

EM12«0B§i2+0i!S12 
E M 2 ls i 2 * ' ' 

EM22»0SS22+0DS22 

RES24»feMii*Y24«Y24*2,#£Ml2«y24#Y42*EM22*Y42*Y42 
P24» 'P<2)#P(4)#<TRA24+TPA42+RES24) 

C FOR P3 A!^0 ^>4‘ 

TRA34«CSlil*0DSll+CSU2«0DS12*DSIll*0CSll*DSI12#0CS12-2,, 

TRA43«CS.U2#0DSl2*CSl22tODS22+DSU24OcSl2+I}SU2#0cS22-2, 

FMii«o&SiitObSii " ■ ■ " ■ ' 

FM12«OCSi2+d0Si2 

FM2laFMi2 

FM22«0CS22+p0S22 

Rgg34*FMii#V34tY34#2.#FMa.2#Y34#Y43*FM22*Y43#Y43 

P34« ~PU)«PT4H<TRA344TRA43+RES34) ' 

Ft ,St (Pi24Pi3^P14+P23*P24i-P34) 

RETURN " 

END , 

JIBFTC CONST 

SUaRbUTiNl C0NST{N,M,K,X,a,H,I) 

DlHiN$|bN X(i>i4)#6(4),H(4) 

Ha)»3p,0 

QC§)« 0 t 

H(2)»30,0 

G«3)»0,, 

H(3)«30,0 

Gt4)s0, 

H(4)s90,0 

X(I»4>sXU#l)'*X(I»2)tX(I,3) 

RETURN ■ ■ ■ 

END 

ilBFTC POPRO 

SUBROUTiNt POPROtXf ieV2*IUIF) 

DIMeNSjbN''Xt6,4),Y(4,2j ' ' 

C0HM0N’THEtA14,4),P(4),0SU4|4)»DETl,DET2»DET3,DET4 i 

1 A(3,3),B(2,2),C(3,3)iD{2,2)iE"(2,2r»AF{2#2hAGt4,4) 

2 AH(2r2} 

integer ' u 
XA*XtlEVS»l) 

X 8 sX<IEV 2 , 2 ? 

XCaX(IEV2i3) 

U®* 2 ' * 

WRITEU*20 )t (THETAd J»ii4)flsl,4> 

WRlTEUiEor « J0SUr»U)fJ»l#4),Ial',4) 

WRlTe<6,20j D6Tl|DET2,D£f3,DET4 
WRnBt6|20) (P(t ),^»l,4r 
20 F0RMAT(2X»’B(aJ(*El2iS))" 

C MODEi. I ' A«BaC 

Ya. is -T HE TA a ,1 )* XA 

Yl3«THETAU'#2)#x8*TH£tA(l,3)#XC 

Y12»-Yii-Yi3' 

iniDIF;Lfc,5) UC TC 35 

YtlfUSYli^O, 6321*0, 003 

Y(li 2) aYlitU ,6 534«0, OQ'i^ 



Y{1*2)s riH£TA(3,lV»XA»^TMTA(3|2)wxB) +0, 6542»0 .0 03 
Y133= CTBETAOf 2)«XB) +0, 9832«0 ,003 
36 OlFll*Y(lfl)-YU 
Diri2=Y(i,2)*Yi2 
CONl«l,0/'t(!)ETl*#0,5)#(6|28##y)) 

pospi*cbfii*^XP<Posi) 

3 HODEl 2 A'nBwO' ' 

y21s-THETS(ii»i)*XAtTHETAU#2)#XB 
Y23sTHlTAr2,3)#XB ' ' 

Y22*-Y21-t23' ' 

Y<2,2}BY<ii2) 

DIF21*Y{2a)-Y21 

DlF22»Yi2i2j-V22 

CON2sl,0/{(C5t2*tO,5}#{6,28»fU)) . 

POS2«-O,5ii0S|<2,l>#l!tF21#Djr2lt2,#OSM2|2)#i)IF23.#DlF22*OSU2,4)*» 

1 IF22*Dir22) . • 

PO$P2«CbN2«EXP(POS2) 

: MO0et. 3 ■ a-Bp-c' 

Y31»«THBTA(3»i)#XA 
Y33PTHeTAt3f2MXB 
Y32#«Y3i»t33 

YC3,8)n'«i#a^ 

DlF3J«f’(3iii-Y3i 

DlF32«Yi3#2j-Y3S 

CON3s3;',b?( (0gtl##Oi5)tt6,28»*U)) . 

POS3s-0,5#^eSI(3il)tPlF3i#ClF3l*2,fbSI<3|2)#OIF3l#DlF32*OSn3|4}»D 

I |F32#DII^‘3^> ' ' ' ■' ' ■ * •■ 

POSP3s0dN3*6XPCPCS3) 

: M0DEl“ 4' 'A=lii=c ■ 

Y4l8-THeTA(4»i>#XA*THETA<4»2)«XB 

Y43*THeTA(4,3j#XB<-tHeTA(4,4j*XC 

Y42*-Y41-Y43" ' ■ 

Y<4#l)»YCi,i) 

Y(4»2)sYU#2j 
DIF41aY'(4;i)»Y41 
UIF42sYi4,2j».V42 , 

CON4sl,O?((0Et4#fO,5)«{6,28#*U)) , 

POS49«»Q,&«iOSr{4ia>*CIF41#0IF4l*2,«OSH4i2)#Bir41#DIF42*OSI(4#4)#D 

llF42«DlF42j' * ' ' ” ' 

POSP4aO0fj44EXP(POS4) 

D£NOHspUT4FbSPi+PC2)«POSP2*P<3)#POSP3^P{4)#POSP4 

PR0BlaP(i?4F&SPl/fiEN0M ' 

PR0B2«P(2)'*PbSF2/pEfv0M 
PR0B3sPi3)#P0SPS/D&N0M . 

PROB4sPUJ«Pb§f»4/0lNdH 
WRITE(6,ai} ■XAf)(B,XCfy(l,2),Yi33 
22 F0RMAt(5FK3)' ■■'•■■• 

21 F0RMAT(S(6xiEl2#5)> 

WRITEUi2iOrFRCBl#PROP2,PROB3,PROH , 

210 FORMAtU5X»4PftOBi*‘4iFa*0»*PRO62«4,FI,S##PROB3»#,F8#5i#PROB4»#,F8f5 

1 ) ■ ' ” ■ ' , 

RETURN 

END 

lENTRY 



